Skip to main content

Section 1.13 SciPy's IVP solvers

Rather than implementing by hand a IVP method to solve a system of \(n\) ODEs, one can choose to use SciPy's own IVP solver. The solver is named scipy.integrate.solve_ivp and takes in input, in its minimal version, three arguments:
  1. A \(n\)-components vector function, called ode in the code below, that returns the righ hand side of the first-order ODE system (if the problem has higher order, the order must e brought to 1 with the usual trick);
  2. A 2-component array (2-tuple) containg the initial and final time;
  3. A \(n\)-components array containing the initial conditions.

For instance, the line

solve_ivp(lambda t, x: t-x, [0, 15], [2])

evaluates a numerical approximation of the IVP
\begin{equation*} \dot x=t-x,\;t_0=0, \;t_f=15,\;x(t_0)=2. \end{equation*}
Important optional inputs are:
  1. \(\hbox{The ODE algorithm to be used. Possible choices:}\)
    'RK45'
    an explicit Runge-Kutta variable-stepsize algorithm of global order 5 where the stepsize is detitleined by running at every step both RK4 and RK5.
    'RK23'
    an explicit Runge-Kutta variable-stepsize algorithm of global order 3 where the stepsize is detitleined by running at every step both RK2 and RK3.
    'DOP853'
    an explicit Runge-Kutta variable-stepsize algorithm of global order 8.
    'Radau'
    an implicit Runge-Kutta variable-stepsize algorithm of global order 5.
    'BFD'
    an implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation.
    'LSODA'
    a multistep method with automatic stiffness detection and switching.
  2. A relative tolerance rtol and an absolute tolerance atol. The solver will keep the local error estimates below atol+rtol*abs(x).
  3. An array t_eval containing particular times at which the user wants to have an evaluation of the solution.
  4. An option dense_output. It must be set to True or it will be impossible to plot a smooth version of the solution (see the code below).
For instance, the line

solve_ivp(lambda t, x: t-x, [0, 15], [2], method = 'DOP853', t_eval=[3.3 11.5], rtol=1e-6)

evaluates a numerical solution to the same IVP above using the DOP853 algorihtm, making sure to obtain a relative error not larger than \(10^{-6}\) at each step and evaluating, in particular, \(x(3.3)\) and \(x(11.5)\text{.}\)

Below we solve our good old IVP problem
\begin{equation*} \dot x=-x\cos t,\;x(0)=1,\;t_f=25. \end{equation*}

And now we solve our 2-dimensional physical pendulum problem
\begin{equation*} \begin{cases} \dot x(t) = y(t),\cr \dot y(t) = -\sin(x(t))\cr \end{cases} \end{equation*}
with initial conditions \(x(0)=1, y(0)=0\) and \(t_f=16.5\text{.}\)
Other Python's ODE Solvers. There is a plethora of other ODE solvers for Python. You can find a thorough list with examples at this url: Ordinary differential equation solvers in Python.