scipy.integrate.solve_ivp¶

scipy.integrate.
solve_ivp
(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)[source]¶ Solve an initial value problem for a system of ODEs.
This function numerically integrates a system of ordinary differential equations given an initial value:
dy / dt = f(t, y) y(t0) = y0
Here t is a onedimensional independent variable (time), y(t) is an ndimensional vectorvalued function (state), and an ndimensional vectorvalued function f(t, y) determines the differential equations. The goal is to find y(t) approximately satisfying the differential equations, given an initial value y(t0)=y0.
Some of the solvers support integration in the complex domain, but note that for stiff ODE solvers, the righthand side must be complexdifferentiable (satisfy CauchyRiemann equations [11]). To solve a problem in the complex domain, pass y0 with a complex data type. Another option is always to rewrite your problem for real and imaginary parts separately.
Parameters:  fun : callable
Righthand side of the system. The calling signature is
fun(t, y)
. Heret
is a scalar, and there are two options for the ndarrayy
: It can either have shape (n,); thenfun
must return array_like with shape (n,). Alternatively it can have shape (n, k); thenfun
must return an array_like with shape (n, k), i.e. each column corresponds to a single column iny
. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows a faster approximation of the Jacobian by finite differences (required for stiff solvers). t_span : 2tuple of floats
Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf.
 y0 : array_like, shape (n,)
Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial guess is purely real).
 method : string or
OdeSolver
, optional Integration method to use:
 ‘RK45’ (default): Explicit RungeKutta method of order 5(4) [1]. The error is controlled assuming accuracy of the fourthorder method, but steps are taken using the fifthorder accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [2]. Can be applied in the complex domain.
 ‘RK23’: Explicit RungeKutta method of order 3(2) [3]. The error is controlled assuming accuracy of the secondorder method, but steps are taken using the thirdorder accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.
 ‘Radau’: Implicit RungeKutta method of the Radau IIA family of order 5 [4]. The error is controlled with a thirdorder accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output.
 ‘BDF’: Implicit multistep variableorder (1 to 5) method based on a backward differentiation formula for the derivative approximation [5]. The implementation follows the one described in [6]. A quasiconstant step scheme is used and accuracy is enhanced using the NDF modification. Can be applied in the complex domain.
 ‘LSODA’: Adams/BDF method with automatic stiffness detection and switching [7], [8]. This is a wrapper of the Fortran solver from ODEPACK.
You should use the ‘RK45’ or ‘RK23’ method for nonstiff problems and ‘Radau’ or ‘BDF’ for stiff problems [9]. If not sure, first try to run ‘RK45’. If needs unusually many iterations, diverges, or fails, your problem is likely to be stiff and you should use ‘Radau’ or ‘BDF’. ‘LSODA’ can also be a good universal choice, but it might be somewhat less convenient to work with as it wraps old Fortran code.
You can also pass an arbitrary class derived from
OdeSolver
which implements the solver. dense_output : bool, optional
Whether to compute a continuous solution. Default is False.
 t_eval : array_like or None, optional
Times at which to store the computed solution, must be sorted and lie within t_span. If None (default), use points selected by the solver.
 events : callable, list of callables or None, optional
Types of events to track. Each is defined by a continuous function of time and state that becomes zero value in case of an event. Each function must have the signature
event(t, y)
and return a float. The solver will find an accurate value oft
at whichevent(t, y(t)) = 0
using a rootfinding algorithm. Additionally eachevent
function might have the following attributes: terminal: bool, whether to terminate integration if this event occurs. Implicitly False if not assigned.
 direction: float, direction of a zero crossing. If direction is positive, event must go from negative to positive, and vice versa if direction is negative. If 0, then either direction will count. Implicitly 0 if not assigned.
You can assign attributes like
event.terminal = True
to any function in Python. If None (default), events won’t be tracked. vectorized : bool, optional
Whether fun is implemented in a vectorized fashion. Default is False.
 options
Options passed to a chosen solver. All options available for already implemented solvers are listed below.
 max_step : float, optional
Maximum allowed step size. Default is np.inf, i.e. the step size is not bounded and determined solely by the solver.
 rtol, atol : float and array_like, optional
Relative and absolute tolerances. The solver keeps the local error estimates less than
atol + rtol * abs(y)
. Here rtol controls a relative accuracy (number of correct digits). But if a component of y is approximately below atol, the error only needs to fall within the same atol threshold, and the number of correct digits is not guaranteed. If components of y have different scales, it might be beneficial to set different atol values for different components by passing array_like with shape (n,) for atol. Default values are 1e3 for rtol and 1e6 for atol. jac : {None, array_like, sparse_matrix, callable}, optional
Jacobian matrix of the righthand side of the system with respect to y, required by the ‘Radau’, ‘BDF’ and ‘LSODA’ method. The Jacobian matrix has shape (n, n) and its element (i, j) is equal to
d f_i / d y_j
. There are three ways to define the Jacobian: If array_like or sparse_matrix, the Jacobian is assumed to be constant. Not supported by ‘LSODA’.
 If callable, the Jacobian is assumed to depend on both
t and y; it will be called as
jac(t, y)
as necessary. For the ‘Radau’ and ‘BDF’ methods, the return value might be a sparse matrix.  If None (default), the Jacobian will be approximated by finite differences.
It is generally recommended to provide the Jacobian rather than relying on a finitedifference approximation.
 jac_sparsity : {None, array_like, sparse matrix}, optional
Defines a sparsity structure of the Jacobian matrix for a finitedifference approximation. Its shape must be (n, n). This argument is ignored if jac is not None. If the Jacobian has only few nonzero elements in each row, providing the sparsity structure will greatly speed up the computations [10]. A zero entry means that a corresponding element in the Jacobian is always zero. If None (default), the Jacobian is assumed to be dense. Not supported by ‘LSODA’, see lband and uband instead.
 lband, uband : int or None
Parameters defining the bandwidth of the Jacobian for the ‘LSODA’ method, i.e.,
jac[i, j] != 0 only for i  lband <= j <= i + uband
. Setting these requires your jac routine to return the Jacobian in the packed format: the returned array must haven
columns anduband + lband + 1
rows in which Jacobian diagonals are written. Specificallyjac_packed[uband + i  j , j] = jac[i, j]
. The same format is used inscipy.linalg.solve_banded
(check for an illustration). These parameters can be also used withjac=None
to reduce the number of Jacobian elements estimated by finite differences. min_step, first_step : float, optional
The minimum allowed step size and the initial step size respectively for ‘LSODA’ method. By default min_step is zero and first_step is selected automatically.
Returns:  Bunch object with the following fields defined:
 t : ndarray, shape (n_points,)
Time points.
 y : ndarray, shape (n, n_points)
Values of the solution at t.
 sol :
OdeSolution
or None Found solution as
OdeSolution
instance; None if dense_output was set to False. t_events : list of ndarray or None
Contains for each event type a list of arrays at which an event of that type event was detected. None if events was None.
 nfev : int
Number of evaluations of the righthand side.
 njev : int
Number of evaluations of the Jacobian.
 nlu : int
Number of LU decompositions.
 status : int
Reason for algorithm termination:
 1: Integration step failed.
 0: The solver successfully reached the end of tspan.
 1: A termination event occurred.
 message : string
Humanreadable description of the termination reason.
 success : bool
True if the solver reached the interval end or a termination event occurred (
status >= 0
).
References
[1] (1, 2) J. R. Dormand, P. J. Prince, “A family of embedded RungeKutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 1926, 1980. [2] (1, 2) L. W. Shampine, “Some Practical RungeKutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135150, 1986. [3] (1, 2) P. Bogacki, L.F. Shampine, “A 3(2) Pair of RungeKutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321325, 1989. [4] (1, 2) E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems”, Sec. IV.8. [5] (1, 2) Backward Differentiation Formula on Wikipedia. [6] (1, 2) L. F. Shampine, M. W. Reichelt, “THE MATLAB ODE SUITE”, SIAM J. SCI. COMPUTE., Vol. 18, No. 1, pp. 122, January 1997. [7] (1, 2) A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 5564, 1983. [8] (1, 2) L. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations”, SIAM Journal on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136148, 1983. [9] (1, 2) Stiff equation on Wikipedia. [10] (1, 2) A. Curtis, M. J. D. Powell, and J. Reid, “On the estimation of sparse Jacobian matrices”, Journal of the Institute of Mathematics and its Applications, 13, pp. 117120, 1974. [11] (1, 2) CauchyRiemann equations on Wikipedia. Examples
Basic exponential decay showing automatically chosen time points.
>>> from scipy.integrate import solve_ivp >>> def exponential_decay(t, y): return 0.5 * y >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8]) >>> print(sol.t) [ 0. 0.11487653 1.26364188 3.06061781 4.85759374 6.65456967 8.4515456 10. ] >>> print(sol.y) [[2. 1.88836035 1.06327177 0.43319312 0.17648948 0.0719045 0.02929499 0.01350938] [4. 3.7767207 2.12654355 0.86638624 0.35297895 0.143809 0.05858998 0.02701876] [8. 7.5534414 4.25308709 1.73277247 0.7059579 0.287618 0.11717996 0.05403753]]
Specifying points where the solution is desired.
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8], ... t_eval=[0, 1, 2, 4, 10]) >>> print(sol.t) [ 0 1 2 4 10] >>> print(sol.y) [[2. 1.21305369 0.73534021 0.27066736 0.01350938] [4. 2.42610739 1.47068043 0.54133472 0.02701876] [8. 4.85221478 2.94136085 1.08266944 0.05403753]]
Cannon fired upward with terminal event upon impact. The
terminal
anddirection
fields of an event are applied by monkey patching a function. Herey[0]
is position andy[1]
is velocity. The projectile starts at position 0 with velocity +10. Note that the integration never reaches t=100 because the event is terminal.>>> def upward_cannon(t, y): return [y[1], 0.5] >>> def hit_ground(t, y): return y[1] >>> hit_ground.terminal = True >>> hit_ground.direction = 1 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground) >>> print(sol.t_events) [array([ 20.])] >>> print(sol.t) [0.00000000e+00 9.99900010e05 1.09989001e03 1.10988901e02 1.11088891e01 1.11098890e+00 1.11099890e+01 2.00000000e+01]