How to solve a single differential equation

12.1 A single ODE

Here is an example of a first order ordinary differential equation (ODE):

$\frac{dy}{dt}=0.9 y^2\cos t$

The ODE is a so-called initial value problem (IVP) whose solution is well-defined once the initial value condition is specified. That may e.g. be that the point:

$(t,y)=(0,1)$

is part of the solution.

Assume we want to solve the ODE for $t$ -values fulfilling $0\le t\le 3$ , we may start by declaring a list (or tuple or NumPy array) with the initial and final $t$ -values:

tinit = 0
tfinal = 3
trange = [tinit,tfinal]

and a list (or NumPy array, but not a tuple) with the $y$ -value dictated by the initial value condition:

yinit = [1]

At this point, it may seem unnecessary to declare this value as a list, but that will turn out logical once we are to solve several differential equations at the same time.

As a final preparation, let us define a function that given $t$ and $y$ evaluates $\displaystyle\frac{dy}{dt}$ :

import numpy as np
dydt = lambda t, y: 0.9 * y**2 * np.cos(t)

We are now ready to import solve_ivp from scipy.integrate and call it like this:

from scipy.integrate import solve_ivp

mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
ys = mysol.y[0]

There are three arguments for solve_ivp:

dydt is the function that implements the differential equation
trange is a list with the starting and ending $t$ -values for which the differential equation should be solved. Importantly, the first element of trange has to conform with that of the initial values.
yinit is a list with the $y$ -value dictated by the initial values.

The ts = mysol.t and ys = mysol.y[0] following the call of solve_ivp select from the output of the function, the $t$ and $y$ values for which the differential equation was actually solved.

These $(t,y)$ values may then be plotted with:

import matplotlib.pyplot as plt

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(ts,ys,'o',linestyle='--')
ax.grid()
ax.set_ylim([0,11])
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')

Figure 1

12.2 ODE: More datapoints

Using solve_ivp to solve a differential equation, the user may specify for which $t$ -values, the solution is sought. That is done by adding a keyword-argument t_eval with a Python list (or NumPy array) containing the $t$ values:

ts = np.linspace(tinit, tfinal, 100)
mysol = solve_ivp(dydt,
                  trange,
                  yinit,
                  t_eval=ts)
ts = mysol.t
ys = mysol.y[0]

ax.plot(ts,ys)
fig

Figure 2

12.3 ODE: Backwards integration

Assume the differential equation is to be solved with this initial value condition:

$(t,y)=(9,\frac{5}{4})$

and for $t$ -values fulfilling $4\le t\le 9$ , i.e. for $t$ values before that of the initial condition. Well, it makes no difference to solve_ivp, that the integration now takes place from large $t$ values to smaller ones. Only, it is important that the first $t$ -value given conforms with that of the initial value condition. Thus, the trange and yinit must be changed into:

tinit = 9
tfinal = 4
trange = (tinit,tfinal)
yinit = [5/4]

After which, the same call of solve_ivp can be made and the result be plotted:

mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
ys = mysol.y[0]

ax.plot(ts,ys,'o',linestyle='--')
ax.set_xlim([-1,10])
fig

Note how trange is defined so that the first value is the one conforming with the new initial values.

Figure 3

12.4 ODE: Accuracy

The differential equation above can be solve by separation of the variables and integration. It turns out to have the exact solution:

$y(t) = \frac{10 y_0}{10-9 y_0(\sin(t)-\sin(t_0)}.$

Using initial values, $(t_0,y_0)=(0,1)$ , and solving until $t=\pi$ one should thus get:

$y(t=\pi) = \frac{10}{10-9\left(\sin(\pi)-\sin(0)\right)} = \frac{10}{10} = 1$

as the final value for $y$ . Lets inspect that:

tinit = 0
tfinal = np.pi
trange = [tinit,tfinal]
yinit = [1]
mysol = solve_ivp(dydt, trange, yinit)

ys = mysol.y[0]

# now print the final value for y, corresponding to y(tfinal) = y(t = pi)
ys[-1]

1.0107395513007404

Well, that was not exactly the case. This is due to rounding errors made by solve_ivp while solving the differential equation.

One may control the accuracy by setting e.g. a maximum step size for the independent variable, here $t$ , during the integration:

mysol = solve_ivp(dydt, trange, yinit, max_step=1e-1)

ys = mysol.y[0]

ys[-1]

1.0000016575704107

Here the maximum step size was set to $0.1$ and the numerical result for $y(\pi)$ came much closer to exact $1$ .

Using np.logspace:

np.logspace(0,-4,5)

array([1.e+00, 1.e-01, 1.e-02, 1.e-03, 1.e-04])

and a list comprehension we may construct a neat list of max_step-values and $y(\pi)$ -values in this way:

res = [(max_step, solve_ivp(dydt,trange,yinit,max_step=max_step).y[0][-1]) \
             for max_step in np.logspace(0,-4,5)]
res

Perhaps not so easy to read. You might prefer to write it like this:

res = []
for max_step in np.logspace(0,-4,5):
  mysol = solve_ivp(dydt, trange, yinit, max_step=max_step)
  ys = mysol.y[0]
  res.append((max_step, ys[-1]))
res

They both give the exact same result, namely:

[(1.0, 1.0107395513007404),
 (0.1, 1.0000016575704107),
 (0.01, 1.0000000000079718),
 (0.001, 0.9999999999999999),
 (0.0001, 0.9999999999999947)]

i.e. practically $1$ with a sufficiently small value for max_step.

Other arguments can be set when calling solve_ivp. These include the relative tolerance, rtol, and the absolute tolerance, atol. You may find more information here.