In scipy, there are also a basic solver for solving the boundary value problems, that is the scipy.integrate.solve_bvp function. The function solves a first order system of ODEs subject to two-point boundary conditions. The function construction are shown below:

CONSTRUCTION:

Let \(F\) be a function object to the function that computes

\(t\) is a one-dimensional independent variable (time), \(S(t)\) is an n-dimensional vector-valued function (state), and the \(F(t, S(t))\) defines the differential equations. \(S0\) be an initial value for \(S\). The function \(F\) must have the form \(dS = F(t, S)\), although the name does not have to be \(F\). The goal is to find the \(S(t)\) approximately satisfying the differential equations, given the initial value \(S(t0)=S0\).

The way we use the solver to solve the differential equation is: $\(solve\_ivp(fun, t\_span, s0, method = 'RK45', t\_eval=None)\)$

where \(fun\) takes in the function in the right-hand side of the system. \(t\_span\) is the interval of integration \((t0, tf)\), where \(t0\) is the start and \(tf\) is the end of the interval. \(s0\) is the initial state. There are a couple of methods that we can choose, the default is ‘RK45’, which is the explicit Runge-Kutta method of order 5(4). There are other methods you can use as well, see the end of this section for more information. \(t\_eval\) takes in the times at which to store the computed solution, and must be sorted and lie within \(t\_span\).

Let’s try one example below.

The above figure shows the corresponding numerical results. As in the previous example, the difference between the result of solve_ivp and the evaluation of the analytical solution by Python is very small in comparison to the value of the function.

EXAMPLE:

Let the state of a system be defined by \(S(t) = \left[\begin{array}{c} x(t) \\y(t) \end{array}\right]\), and let the evolution of the system be defined by the ODE

Use solve_ivp to solve this ODE for the time interval \([0, 10]\) with an initial value of \(S_0 = \left[\begin{array}{c} 1 \\1 \end{array}\right]\). Plot the solution in (\(x(t), y(t)\)).

In practice, some ODEs have bad behavior known as stiffness. Loosely speaking, stiffness refers to systems that can have very sharp changes in derivative. An example of a stiff system is a bouncing ball, which suddenly changes directions when it hits the ground. Depending on the properties of the ODE you are solving and the desired level of accuracy, you might need to use different methods for solve_ivp. There are many methods that you can choose for the method argument in solve_ivp, take a look of the documentation to understand it more. As suggested by the documentation, you should use the ‘RK45’ or ‘RK23’ method for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems. 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.

< 23.4 Numerical Error and Instability | Contents | 23.6 Summary and Problems >