scipy.integrate._ivp.ivp:solve_ivp

Parameters

fun : callable

Right-hand side of the system: the time derivative of the state y at time t. The calling signature is fun(t, y), where t is a scalar and y is an ndarray with len(y) = len(y0). Additional arguments need to be passed if args is used (see documentation of args argument). fun must return an array of the same shape as y. See vectorized for more information.

t_span : 2-member sequence

Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf. Both t0 and tf must be floats or values interpretable by the float conversion function.

y0 : array_like, shape (n,)

Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial value is purely real).

method : string or `OdeSolver`, optional

Integration method to use:

• 'RK45' (default): Explicit Runge-Kutta method of order 5(4) ^[1]. The error is controlled assuming accuracy of the fourth-order method, but steps are taken using the fifth-order 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 Runge-Kutta method of order 3(2) ^[3]. The error is controlled assuming accuracy of the second-order method, but steps are taken using the third-order accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.

• 'DOP853': Explicit Runge-Kutta method of order 8 ^[13]. Python implementation of the "DOP853" algorithm originally written in Fortran ^[14]. A 7-th order interpolation polynomial accurate to 7-th order is used for the dense output. Can be applied in the complex domain.

• 'Radau': Implicit Runge-Kutta method of the Radau IIA family of order 5 ^[4]. The error is controlled with a third-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 formula for the derivative approximation ^[5]. The implementation follows the one described in ^[6]. A quasi-constant 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.

Explicit Runge-Kutta methods ('RK23', 'RK45', 'DOP853') should be used for non-stiff problems and implicit methods ('Radau', 'BDF') for stiff problems ^[9]. Among Runge-Kutta methods, 'DOP853' is recommended for solving with high precision (low values of rtol and atol).

If not sure, first try to run 'RK45'. If it makes 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.

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.

dense_output : bool, optional

Whether to compute a continuous solution. Default is False.

events : callable, or list of callables, optional

Events to track. If None (default), no events will be tracked. Each event occurs at the zeros of a continuous function of time and state. Each function must have the signature event(t, y) where additional argument have to be passed if args is used (see documentation of args argument). Each function must return a float. The solver will find an accurate value of t at which event(t, y(t)) = 0 using a root-finding algorithm. By default, all zeros will be found. The solver looks for a sign change over each step, so if multiple zero crossings occur within one step, events may be missed. Additionally each event function might have the following attributes:

terminal: bool or int, optional

When boolean, whether to terminate integration if this event occurs. When integral, termination occurs after the specified the number of occurrences of this event. Implicitly False if not assigned.

direction: float, optional

Direction of a zero crossing. If direction is positive, event will only trigger when going from negative to positive, and vice versa if direction is negative. If 0, then either direction will trigger event. Implicitly 0 if not assigned.

You can assign attributes like event.terminal = True to any function in Python.

vectorized : bool, optional

Whether fun can be called in a vectorized fashion. Default is False.

If vectorized is False, fun will always be called with y of shape (n,), where n = len(y0).

If vectorized is True, fun may be called with y of shape (n, k), where k is an integer. In this case, fun must behave such that fun(t, y)[:, i] == fun(t, y[:, i]) (i.e. each column of the returned array is the time derivative of the state corresponding with a column of y).

Setting vectorized=True allows for faster finite difference approximation of the Jacobian by methods 'Radau' and 'BDF', but will result in slower execution for other methods and for 'Radau' and 'BDF' in some circumstances (e.g. small len(y0)).

args : tuple, optional

Additional arguments to pass to the user-defined functions. If given, the additional arguments are passed to all user-defined functions. So if, for example, fun has the signature fun(t, y, a, b, c), then jac (if given) and any event functions must have the same signature, and args must be a tuple of length 3.

**options

Options passed to a chosen solver. All options available for already implemented solvers are listed below.

first_step : float or None, optional

Initial step size. Default is None which means that the algorithm should choose.

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 or 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), while atol controls absolute accuracy (number of correct decimal places). To achieve the desired rtol, set atol to be smaller than the smallest value that can be expected from rtol * abs(y) so that rtol dominates the allowable error. If atol is larger than rtol * abs(y) the number of correct digits is not guaranteed. Conversely, to achieve the desired atol set rtol such that rtol * abs(y) is always smaller than atol. 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 : array_like, sparse_matrix, callable or None, optional

Jacobian matrix of the right-hand 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. Additional arguments have to be passed if args is used (see documentation of args argument). For '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 finite-difference approximation.

jac_sparsity : array_like, sparse matrix or None, optional

Defines a sparsity structure of the Jacobian matrix for a finite- difference approximation. Its shape must be (n, n). This argument is ignored if jac is not None. 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 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, optional

Parameters defining the bandwidth of the Jacobian for the 'LSODA' method, i.e., jac[i, j] != 0 only for i - lband <= j <= i + uband. Default is None. 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 : float, optional

The minimum allowed step size for 'LSODA' method. By default min_step is zero.