Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_odepack_py.py: 29%
31 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1# Author: Travis Oliphant
3__all__ = ['odeint']
5import numpy as np
6from . import _odepack
7from copy import copy
8import warnings
11class ODEintWarning(Warning):
12 pass
15_msgs = {2: "Integration successful.",
16 1: "Nothing was done; the integration time was 0.",
17 -1: "Excess work done on this call (perhaps wrong Dfun type).",
18 -2: "Excess accuracy requested (tolerances too small).",
19 -3: "Illegal input detected (internal error).",
20 -4: "Repeated error test failures (internal error).",
21 -5: "Repeated convergence failures (perhaps bad Jacobian or tolerances).",
22 -6: "Error weight became zero during problem.",
23 -7: "Internal workspace insufficient to finish (internal error).",
24 -8: "Run terminated (internal error)."
25 }
28def odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0,
29 ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0,
30 hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
31 mxords=5, printmessg=0, tfirst=False):
32 """
33 Integrate a system of ordinary differential equations.
35 .. note:: For new code, use `scipy.integrate.solve_ivp` to solve a
36 differential equation.
38 Solve a system of ordinary differential equations using lsoda from the
39 FORTRAN library odepack.
41 Solves the initial value problem for stiff or non-stiff systems
42 of first order ode-s::
44 dy/dt = func(y, t, ...) [or func(t, y, ...)]
46 where y can be a vector.
48 .. note:: By default, the required order of the first two arguments of
49 `func` are in the opposite order of the arguments in the system
50 definition function used by the `scipy.integrate.ode` class and
51 the function `scipy.integrate.solve_ivp`. To use a function with
52 the signature ``func(t, y, ...)``, the argument `tfirst` must be
53 set to ``True``.
55 Parameters
56 ----------
57 func : callable(y, t, ...) or callable(t, y, ...)
58 Computes the derivative of y at t.
59 If the signature is ``callable(t, y, ...)``, then the argument
60 `tfirst` must be set ``True``.
61 y0 : array
62 Initial condition on y (can be a vector).
63 t : array
64 A sequence of time points for which to solve for y. The initial
65 value point should be the first element of this sequence.
66 This sequence must be monotonically increasing or monotonically
67 decreasing; repeated values are allowed.
68 args : tuple, optional
69 Extra arguments to pass to function.
70 Dfun : callable(y, t, ...) or callable(t, y, ...)
71 Gradient (Jacobian) of `func`.
72 If the signature is ``callable(t, y, ...)``, then the argument
73 `tfirst` must be set ``True``.
74 col_deriv : bool, optional
75 True if `Dfun` defines derivatives down columns (faster),
76 otherwise `Dfun` should define derivatives across rows.
77 full_output : bool, optional
78 True if to return a dictionary of optional outputs as the second output
79 printmessg : bool, optional
80 Whether to print the convergence message
81 tfirst : bool, optional
82 If True, the first two arguments of `func` (and `Dfun`, if given)
83 must ``t, y`` instead of the default ``y, t``.
85 .. versionadded:: 1.1.0
87 Returns
88 -------
89 y : array, shape (len(t), len(y0))
90 Array containing the value of y for each desired time in t,
91 with the initial value `y0` in the first row.
92 infodict : dict, only returned if full_output == True
93 Dictionary containing additional output information
95 ======= ============================================================
96 key meaning
97 ======= ============================================================
98 'hu' vector of step sizes successfully used for each time step
99 'tcur' vector with the value of t reached for each time step
100 (will always be at least as large as the input times)
101 'tolsf' vector of tolerance scale factors, greater than 1.0,
102 computed when a request for too much accuracy was detected
103 'tsw' value of t at the time of the last method switch
104 (given for each time step)
105 'nst' cumulative number of time steps
106 'nfe' cumulative number of function evaluations for each time step
107 'nje' cumulative number of jacobian evaluations for each time step
108 'nqu' a vector of method orders for each successful step
109 'imxer' index of the component of largest magnitude in the
110 weighted local error vector (e / ewt) on an error return, -1
111 otherwise
112 'lenrw' the length of the double work array required
113 'leniw' the length of integer work array required
114 'mused' a vector of method indicators for each successful time step:
115 1: adams (nonstiff), 2: bdf (stiff)
116 ======= ============================================================
118 Other Parameters
119 ----------------
120 ml, mu : int, optional
121 If either of these are not None or non-negative, then the
122 Jacobian is assumed to be banded. These give the number of
123 lower and upper non-zero diagonals in this banded matrix.
124 For the banded case, `Dfun` should return a matrix whose
125 rows contain the non-zero bands (starting with the lowest diagonal).
126 Thus, the return matrix `jac` from `Dfun` should have shape
127 ``(ml + mu + 1, len(y0))`` when ``ml >=0`` or ``mu >=0``.
128 The data in `jac` must be stored such that ``jac[i - j + mu, j]``
129 holds the derivative of the `i`th equation with respect to the `j`th
130 state variable. If `col_deriv` is True, the transpose of this
131 `jac` must be returned.
132 rtol, atol : float, optional
133 The input parameters `rtol` and `atol` determine the error
134 control performed by the solver. The solver will control the
135 vector, e, of estimated local errors in y, according to an
136 inequality of the form ``max-norm of (e / ewt) <= 1``,
137 where ewt is a vector of positive error weights computed as
138 ``ewt = rtol * abs(y) + atol``.
139 rtol and atol can be either vectors the same length as y or scalars.
140 Defaults to 1.49012e-8.
141 tcrit : ndarray, optional
142 Vector of critical points (e.g., singularities) where integration
143 care should be taken.
144 h0 : float, (0: solver-determined), optional
145 The step size to be attempted on the first step.
146 hmax : float, (0: solver-determined), optional
147 The maximum absolute step size allowed.
148 hmin : float, (0: solver-determined), optional
149 The minimum absolute step size allowed.
150 ixpr : bool, optional
151 Whether to generate extra printing at method switches.
152 mxstep : int, (0: solver-determined), optional
153 Maximum number of (internally defined) steps allowed for each
154 integration point in t.
155 mxhnil : int, (0: solver-determined), optional
156 Maximum number of messages printed.
157 mxordn : int, (0: solver-determined), optional
158 Maximum order to be allowed for the non-stiff (Adams) method.
159 mxords : int, (0: solver-determined), optional
160 Maximum order to be allowed for the stiff (BDF) method.
162 See Also
163 --------
164 solve_ivp : solve an initial value problem for a system of ODEs
165 ode : a more object-oriented integrator based on VODE
166 quad : for finding the area under a curve
168 Examples
169 --------
170 The second order differential equation for the angle `theta` of a
171 pendulum acted on by gravity with friction can be written::
173 theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0
175 where `b` and `c` are positive constants, and a prime (') denotes a
176 derivative. To solve this equation with `odeint`, we must first convert
177 it to a system of first order equations. By defining the angular
178 velocity ``omega(t) = theta'(t)``, we obtain the system::
180 theta'(t) = omega(t)
181 omega'(t) = -b*omega(t) - c*sin(theta(t))
183 Let `y` be the vector [`theta`, `omega`]. We implement this system
184 in Python as:
186 >>> import numpy as np
187 >>> def pend(y, t, b, c):
188 ... theta, omega = y
189 ... dydt = [omega, -b*omega - c*np.sin(theta)]
190 ... return dydt
191 ...
193 We assume the constants are `b` = 0.25 and `c` = 5.0:
195 >>> b = 0.25
196 >>> c = 5.0
198 For initial conditions, we assume the pendulum is nearly vertical
199 with `theta(0)` = `pi` - 0.1, and is initially at rest, so
200 `omega(0)` = 0. Then the vector of initial conditions is
202 >>> y0 = [np.pi - 0.1, 0.0]
204 We will generate a solution at 101 evenly spaced samples in the interval
205 0 <= `t` <= 10. So our array of times is:
207 >>> t = np.linspace(0, 10, 101)
209 Call `odeint` to generate the solution. To pass the parameters
210 `b` and `c` to `pend`, we give them to `odeint` using the `args`
211 argument.
213 >>> from scipy.integrate import odeint
214 >>> sol = odeint(pend, y0, t, args=(b, c))
216 The solution is an array with shape (101, 2). The first column
217 is `theta(t)`, and the second is `omega(t)`. The following code
218 plots both components.
220 >>> import matplotlib.pyplot as plt
221 >>> plt.plot(t, sol[:, 0], 'b', label='theta(t)')
222 >>> plt.plot(t, sol[:, 1], 'g', label='omega(t)')
223 >>> plt.legend(loc='best')
224 >>> plt.xlabel('t')
225 >>> plt.grid()
226 >>> plt.show()
227 """
229 if ml is None:
230 ml = -1 # changed to zero inside function call
231 if mu is None:
232 mu = -1 # changed to zero inside function call
234 dt = np.diff(t)
235 if not ((dt >= 0).all() or (dt <= 0).all()):
236 raise ValueError("The values in t must be monotonically increasing "
237 "or monotonically decreasing; repeated values are "
238 "allowed.")
240 t = copy(t)
241 y0 = copy(y0)
242 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
243 full_output, rtol, atol, tcrit, h0, hmax, hmin,
244 ixpr, mxstep, mxhnil, mxordn, mxords,
245 int(bool(tfirst)))
246 if output[-1] < 0:
247 warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."
248 warnings.warn(warning_msg, ODEintWarning)
249 elif printmessg:
250 warning_msg = _msgs[output[-1]]
251 warnings.warn(warning_msg, ODEintWarning)
253 if full_output:
254 output[1]['message'] = _msgs[output[-1]]
256 output = output[:-1]
257 if len(output) == 1:
258 return output[0]
259 else:
260 return output