Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_ivp/lsoda.py: 19%
57 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
2from scipy.integrate import ode
3from .common import validate_tol, validate_first_step, warn_extraneous
4from .base import OdeSolver, DenseOutput
7class LSODA(OdeSolver):
8 """Adams/BDF method with automatic stiffness detection and switching.
10 This is a wrapper to the Fortran solver from ODEPACK [1]_. It switches
11 automatically between the nonstiff Adams method and the stiff BDF method.
12 The method was originally detailed in [2]_.
14 Parameters
15 ----------
16 fun : callable
17 Right-hand side of the system: the time derivative of the state ``y``
18 at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
19 scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must
20 return an array of the same shape as ``y``. See `vectorized` for more
21 information.
22 t0 : float
23 Initial time.
24 y0 : array_like, shape (n,)
25 Initial state.
26 t_bound : float
27 Boundary time - the integration won't continue beyond it. It also
28 determines the direction of the integration.
29 first_step : float or None, optional
30 Initial step size. Default is ``None`` which means that the algorithm
31 should choose.
32 min_step : float, optional
33 Minimum allowed step size. Default is 0.0, i.e., the step size is not
34 bounded and determined solely by the solver.
35 max_step : float, optional
36 Maximum allowed step size. Default is np.inf, i.e., the step size is not
37 bounded and determined solely by the solver.
38 rtol, atol : float and array_like, optional
39 Relative and absolute tolerances. The solver keeps the local error
40 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
41 relative accuracy (number of correct digits), while `atol` controls
42 absolute accuracy (number of correct decimal places). To achieve the
43 desired `rtol`, set `atol` to be smaller than the smallest value that
44 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
45 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
46 number of correct digits is not guaranteed. Conversely, to achieve the
47 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
48 than `atol`. If components of y have different scales, it might be
49 beneficial to set different `atol` values for different components by
50 passing array_like with shape (n,) for `atol`. Default values are
51 1e-3 for `rtol` and 1e-6 for `atol`.
52 jac : None or callable, optional
53 Jacobian matrix of the right-hand side of the system with respect to
54 ``y``. The Jacobian matrix has shape (n, n) and its element (i, j) is
55 equal to ``d f_i / d y_j``. The function will be called as
56 ``jac(t, y)``. If None (default), the Jacobian will be
57 approximated by finite differences. It is generally recommended to
58 provide the Jacobian rather than relying on a finite-difference
59 approximation.
60 lband, uband : int or None
61 Parameters defining the bandwidth of the Jacobian,
62 i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``. Setting
63 these requires your jac routine to return the Jacobian in the packed format:
64 the returned array must have ``n`` columns and ``uband + lband + 1``
65 rows in which Jacobian diagonals are written. Specifically
66 ``jac_packed[uband + i - j , j] = jac[i, j]``. The same format is used
67 in `scipy.linalg.solve_banded` (check for an illustration).
68 These parameters can be also used with ``jac=None`` to reduce the
69 number of Jacobian elements estimated by finite differences.
70 vectorized : bool, optional
71 Whether `fun` may be called in a vectorized fashion. False (default)
72 is recommended for this solver.
74 If ``vectorized`` is False, `fun` will always be called with ``y`` of
75 shape ``(n,)``, where ``n = len(y0)``.
77 If ``vectorized`` is True, `fun` may be called with ``y`` of shape
78 ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave
79 such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
80 the returned array is the time derivative of the state corresponding
81 with a column of ``y``).
83 Setting ``vectorized=True`` allows for faster finite difference
84 approximation of the Jacobian by methods 'Radau' and 'BDF', but
85 will result in slower execution for this solver.
87 Attributes
88 ----------
89 n : int
90 Number of equations.
91 status : string
92 Current status of the solver: 'running', 'finished' or 'failed'.
93 t_bound : float
94 Boundary time.
95 direction : float
96 Integration direction: +1 or -1.
97 t : float
98 Current time.
99 y : ndarray
100 Current state.
101 t_old : float
102 Previous time. None if no steps were made yet.
103 nfev : int
104 Number of evaluations of the right-hand side.
105 njev : int
106 Number of evaluations of the Jacobian.
108 References
109 ----------
110 .. [1] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
111 Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
112 pp. 55-64, 1983.
113 .. [2] L. Petzold, "Automatic selection of methods for solving stiff and
114 nonstiff systems of ordinary differential equations", SIAM Journal
115 on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
116 1983.
117 """
118 def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,
119 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,
120 uband=None, vectorized=False, **extraneous):
121 warn_extraneous(extraneous)
122 super().__init__(fun, t0, y0, t_bound, vectorized)
124 if first_step is None:
125 first_step = 0 # LSODA value for automatic selection.
126 else:
127 first_step = validate_first_step(first_step, t0, t_bound)
129 first_step *= self.direction
131 if max_step == np.inf:
132 max_step = 0 # LSODA value for infinity.
133 elif max_step <= 0:
134 raise ValueError("`max_step` must be positive.")
136 if min_step < 0:
137 raise ValueError("`min_step` must be nonnegative.")
139 rtol, atol = validate_tol(rtol, atol, self.n)
141 solver = ode(self.fun, jac)
142 solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,
143 min_step=min_step, first_step=first_step,
144 lband=lband, uband=uband)
145 solver.set_initial_value(y0, t0)
147 # Inject t_bound into rwork array as needed for itask=5.
148 solver._integrator.rwork[0] = self.t_bound
149 solver._integrator.call_args[4] = solver._integrator.rwork
151 self._lsoda_solver = solver
153 def _step_impl(self):
154 solver = self._lsoda_solver
155 integrator = solver._integrator
157 # From lsoda.step and lsoda.integrate itask=5 means take a single
158 # step and do not go past t_bound.
159 itask = integrator.call_args[2]
160 integrator.call_args[2] = 5
161 solver._y, solver.t = integrator.run(
162 solver.f, solver.jac or (lambda: None), solver._y, solver.t,
163 self.t_bound, solver.f_params, solver.jac_params)
164 integrator.call_args[2] = itask
166 if solver.successful():
167 self.t = solver.t
168 self.y = solver._y
169 # From LSODA Fortran source njev is equal to nlu.
170 self.njev = integrator.iwork[12]
171 self.nlu = integrator.iwork[12]
172 return True, None
173 else:
174 return False, 'Unexpected istate in LSODA.'
176 def _dense_output_impl(self):
177 iwork = self._lsoda_solver._integrator.iwork
178 rwork = self._lsoda_solver._integrator.rwork
180 order = iwork[14]
181 h = rwork[11]
182 yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
183 (self.n, order + 1), order='F').copy()
185 return LsodaDenseOutput(self.t_old, self.t, h, order, yh)
188class LsodaDenseOutput(DenseOutput):
189 def __init__(self, t_old, t, h, order, yh):
190 super().__init__(t_old, t)
191 self.h = h
192 self.yh = yh
193 self.p = np.arange(order + 1)
195 def _call_impl(self, t):
196 if t.ndim == 0:
197 x = ((t - self.t) / self.h) ** self.p
198 else:
199 x = ((t - self.t) / self.h) ** self.p[:, None]
201 return np.dot(self.yh, x)