Previous topic

scipy.integrate.romb

Next topic

scipy.integrate.RK23

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 1-dimensional independent variable (time), y(t) is an n-dimensional vector-valued function (state) and an n-dimensional vector-valued 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 a complex domain, but note that for stiff ODE solvers the right hand side must be complex differentiable (satisfy Cauchy-Riemann equations [11]). To solve a problem in a complex domain, pass y0 with a complex data type. Another option always available is to rewrite your problem for real and imaginary parts separately.

Parameters:

Parameters:	fun : callable Right-hand side of the system. The calling signature is `fun(t, y)`. Here `t` is a scalar and there are two options for ndarray `y`. It can either have shape (n,), then `fun` must return array_like with shape (n,). Or alternatively it can have shape (n, k), then `fun` must return array_like with shape (n, k), i.e. each column corresponds to a single column in `y`. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows faster approximation of the Jacobian by finite differences (required for stiff solvers). t_span : 2-tuple 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 a 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 Runge-Kutta method of order 5(4) [R68]. The error is controlled assuming 4th order accuracy, but steps are taken using a 5th oder accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [R69]. Can be applied in a complex domain. ‘RK23’: Explicit Runge-Kutta method of order 3(2) [R70]. The error is controlled assuming 2nd order accuracy, but steps are taken using a 3rd oder accurate formula (local extrapolation is done). A cubic Hermit polynomial is used for the dense output. Can be applied in a complex domain. ‘Radau’: Implicit Runge-Kutta method of Radau IIA family of order 5 [R71]. The error is controlled for a 3rd order accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output. ‘BDF’: Implicit multi-step variable order (1 to 5) method based on a Backward Differentiation Formulas for the derivative approximation [R72]. An implementation approach follows the one described in [R73]. A quasi-constant step scheme is used and accuracy enhancement using NDF modification is also implemented. Can be applied in a complex domain. ‘LSODA’: Adams/BDF method with automatic stiffness detection and switching [R74], [R75]. This is a wrapper of the Fortran solver from ODEPACK. You should use ‘RK45’ or ‘RK23’ methods for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems [R76]. If not sure, first try to run ‘RK45’ and if it does unusual many iterations or diverges then 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 an 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 a solver. events : callable, list of callables or None, optional Events to track. Events are defined by functions which take a zero value at a point of an event. Each function must have a signature `event(t, y)` and return float, the solver will find an accurate value of `t` at which `event(t, y(t)) = 0` using a root finding algorithm. Additionally each `event` function might have attributes: terminal: bool, whether to terminate integration if this event occurs. Implicitly False if not assigned. direction: float, direction of crossing a zero. If direction is positive then event must go from negative to positive, and vice-versa if direction is negative. If 0, then either way 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 constructor. All options available for already implemented solvers are listed below. max_step : float, optional Maximum allowed step size. Default is np.inf, i.e. step 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 then 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 1e-3 for rtol and 1e-6 for atol. jac : {None, array_like, sparse_matrix, callable}, optional Jacobian matrix of the right-hand side of the system with respect to y, required by ‘Radau’, ‘BDF’ and ‘LSODA’ methods. The Jacobian matrix has shape (n, n) and its element (i, j) is equal to `d f_i / d y_j`. There are 3 ways to define the Jacobian: If array_like or sparse_matrix, then the Jacobian is assumed to be constant. Not supported by ‘LSODA’. If callable, then the Jacobian is assumed to depend on both t and y, and will be called as `jac(t, y)` as necessary. For ‘Radau’ and ‘BDF’ methods the return value might be a sparse matrix. If None (default), then the Jacobian will be approximated by finite differences. It is generally recommended to provide the Jacobian rather than relying on a finite difference approximation. jac_sparsity : {None, array_like, sparse matrix}, optional Defines a sparsity structure of the Jacobian matrix for a finite difference approximation, its shape must be (n, n). If the Jacobian has only few non-zero 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 identically 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 Jacobian matrix bandwidth for ‘LSODA’ method. The Jacobian bandwidth means that `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 have `n` columns and `uband + lband + 1` rows in which Jacobian diagonals are written. Specifically `jac_packed[uband + i - j , j] = jac[i, j]`. The same format is used in `scipy.linalg.solve_banded` (check for an illustration). These parameters can be also used with `jac=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) Solution values 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 arrays with times at each a corresponding event was detected, the length of the list equals to the number of events. None if events was None. nfev : int Number of the system rhs evaluations. njev : int Number of the Jacobian evaluations. nlu : int Number of LU decompositions. status : int Reason for algorithm termination: -1: Integration step failed. 0: The solver successfully reached the interval end. 1: A termination event occurred. message : string Verbal description of the termination reason. success : bool True if the solver reached the interval end or a termination event occurred (`status >= 0`).

fun : callable

Right-hand side of the system. The calling signature is fun(t, y). Here t is a scalar and there are two options for ndarray y. It can either have shape (n,), then fun must return array_like with shape (n,). Or alternatively it can have shape (n, k), then fun must return array_like with shape (n, k), i.e. each column corresponds to a single column in y. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows faster approximation of the Jacobian by finite differences (required for stiff solvers).

t_span : 2-tuple 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 a 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 Runge-Kutta method of order 5(4) [R68]. The error is controlled assuming 4th order accuracy, but steps are taken using a 5th oder accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [R69]. Can be applied in a complex domain.

‘RK23’: Explicit Runge-Kutta method of order 3(2) [R70]. The error is controlled assuming 2nd order accuracy, but steps are taken using a 3rd oder accurate formula (local extrapolation is done). A cubic Hermit polynomial is used for the dense output. Can be applied in a complex domain.

‘Radau’: Implicit Runge-Kutta method of Radau IIA family of order 5 [R71]. The error is controlled for a 3rd order accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output.

‘BDF’: Implicit multi-step variable order (1 to 5) method based on a Backward Differentiation Formulas for the derivative approximation [R72]. An implementation approach follows the one described in [R73]. A quasi-constant step scheme is used and accuracy enhancement using NDF modification is also implemented. Can be applied in a complex domain.

‘LSODA’: Adams/BDF method with automatic stiffness detection and switching [R74], [R75]. This is a wrapper of the Fortran solver from ODEPACK.

You should use ‘RK45’ or ‘RK23’ methods for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems [R76]. If not sure, first try to run ‘RK45’ and if it does unusual many iterations or diverges then 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 an 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 a solver.

events : callable, list of callables or None, optional

Events to track. Events are defined by functions which take a zero value at a point of an event. Each function must have a signature event(t, y) and return float, the solver will find an accurate value of t at which event(t, y(t)) = 0 using a root finding algorithm. Additionally each event function might have attributes:

terminal: bool, whether to terminate integration if this event occurs. Implicitly False if not assigned.

direction: float, direction of crossing a zero. If direction is positive then event must go from negative to positive, and vice-versa if direction is negative. If 0, then either way 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 constructor. All options available for already implemented solvers are listed below.

max_step : float, optional

Maximum allowed step size. Default is np.inf, i.e. step 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 then 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 1e-3 for rtol and 1e-6 for atol.

jac : {None, array_like, sparse_matrix, callable}, optional

Jacobian matrix of the right-hand side of the system with respect to y, required by ‘Radau’, ‘BDF’ and ‘LSODA’ methods. The Jacobian matrix has shape (n, n) and its element (i, j) is equal to d f_i / d y_j. There are 3 ways to define the Jacobian:

If array_like or sparse_matrix, then the Jacobian is assumed to be constant. Not supported by ‘LSODA’.

If callable, then the Jacobian is assumed to depend on both t and y, and will be called as jac(t, y) as necessary. For ‘Radau’ and ‘BDF’ methods the return value might be a sparse matrix.

If None (default), then the Jacobian will be approximated by finite differences.

It is generally recommended to provide the Jacobian rather than relying on a finite difference approximation.

jac_sparsity : {None, array_like, sparse matrix}, optional

Defines a sparsity structure of the Jacobian matrix for a finite difference approximation, its shape must be (n, n). If the Jacobian has only few non-zero 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 identically 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 Jacobian matrix bandwidth for ‘LSODA’ method. The Jacobian bandwidth means that 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 have n columns and uband + lband + 1 rows in which Jacobian diagonals are written. Specifically jac_packed[uband + i - j , j] = jac[i, j]. The same format is used in scipy.linalg.solve_banded (check for an illustration). These parameters can be also used with jac=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)

Solution values 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 arrays with times at each a corresponding event was detected, the length of the list equals to the number of events. None if events was None.

nfev : int

Number of the system rhs evaluations.

njev : int

Number of the Jacobian evaluations.

nlu : int

Number of LU decompositions.

status : int

Reason for algorithm termination:

-1: Integration step failed.

0: The solver successfully reached the interval end.

1: A termination event occurred.

message : string

Verbal description of the termination reason.

success : bool

True if the solver reached the interval end or a termination event occurred (status >= 0).

References

[R68]

(1, 2) J. R. Dormand, P. J. Prince, “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.

[R69]

(1, 2) L. W. Shampine, “Some Practical Runge-Kutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.

[R70]

(1, 2) P. Bogacki, L.F. Shampine, “A 3(2) Pair of Runge-Kutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.

[R71]

(1, 2) E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”, Sec. IV.8.

[R72]

(1, 2) Backward Differentiation Formula on Wikipedia.

[R73]

(1, 2) L. F. Shampine, M. W. Reichelt, “THE MATLAB ODE SUITE”, SIAM J. SCI. COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.

[R74]

(1, 2) A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 55-64, 1983.

[R75]

(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. 136-148, 1983.

[R76]

(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. 117-120, 1974.

[11]	(1, 2) Cauchy-Riemann 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 and direction fields of an event are applied by monkey patching a function. Here y[0] is position and y[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.99900010e-05   1.09989001e-03   1.10988901e-02
   1.11088891e-01   1.11098890e+00   1.11099890e+01   2.00000000e+01]