Euler's method and the midpoint method for solving ODEs¶
This notebook covers roughly the same ground as this video.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
Euler's method is the most obvious method for solving an ordinary differential equation of the form $dx/dt=f(t,x)$, but it is significantly less accurate than more sophisticated methods that will be discussed later. As input, we need a time interval $[a,b]$ and an initial value $x_0$ for $x$ at time $t=a$, and also a number $n$ of steps. We divide our time interval into subintervals $[t_i,t_{i+1}]$ of length $h=(b-a)/n$, for $0\leq i\lt n$. For $0\leq i\lt n$ we then define $x_{i+1}=x_i+h\,f(t_i,x_i)$; this is an approximation to the value of the solution at time $t_i$.
def euler(f, a, b, x0, n = 100):
"""
Solve dx/dt = f(t,x) for t from a to b using n steps
Return a pair (ts, xs) where ts is the array of values of t
and xs is the array of corresponding values of x.
"""
# ts is an array of n+1 equally spaced numbers from a to b
# There are n+1 numbers, with n gaps between them, each of size (b-a)/n
ts = np.linspace(a, b, n+1)
# If x0 is just a scalar, the line below is equivalent to xs = np.zeros(n+1)
# We use a fancier version for later cases where x0 is a numpy array of arbitrary shape
xs = np.zeros((n+1,) + np.shape(x0))
xs[0] = x0
h = (b-a)/n
for i in range(n):
xs[i+1] = xs[i] + h * f(ts[i], xs[i])
return ts, xs
We now use Euler's method to solve $dx/dt=x$ for $0\leq t\leq 1$ with $x=1$ when $t=0$. The exact solution is of course $x=e^t$; this is plotted in blue. We have also plotted the result of Euler's method with increasing numbers of steps. With $10$ steps the result is quite noticeably different from the true value. With more steps the error gets smaller, but is still rather poor.
f0 = lambda t, x: x
plt.figure(figsize=(20, 6))
plt.ylim(0, 3)
ts = np.linspace(0, 1, 101)
plt.plot(ts, np.exp(ts), label='exact' )
for k in range(10,80,10):
ts, xs = euler(f0, 0, 1, 1, k)
plt.plot(ts, xs, label=f'{k} steps')
plt.legend()
<matplotlib.legend.Legend at 0x21121cd3350>
The midpoint method is better than Euler's method. Suppose we have found a value for $x$ at time $t$, and we want a value at time $t+h$. We first use Euler's method with a step length of $h/2$ to get an estimated value for $x$ at time $t^*=t+h/2$, namely $x^*=x+f(t,x)h/2$. We then use $f(t^*,x^*)$ as the estimated value of $dx/dt$, so the estimated value of $x$ at $t+h$ is $x+f(t^*,x^*)h$.
def midpoint(f, t0, t1, x0, n = 100):
ts = np.linspace(t0, t1, n+1)
# If x0 is just a scalar, the line below is equivalent to xs = np.zeros(n+1)
# We use a fancier version for later cases where x0 is a numpy array of arbitrary shape
xs = np.zeros((n+1,) + np.shape(x0))
xs[0] = x0
h = (t1-t0)/n
for i in range(n):
tt = ts[i] + h/2
xx = xs[i] + h/2 * f(ts[i], xs[i])
xs[i+1] = xs[i] + h * f(tt, xx)
return ts, xs
We again solve $dx/dt$ with $x=1$ when $t=0$, using the midpoint method. We find that even with only $3$ steps, the error is quite small, and it becomes invisible with 12 steps.
plt.figure(figsize=(20, 6))
plt.ylim(0, 3)
ts = np.linspace(0, 1, 101)
plt.plot(ts, np.exp(ts), label='exact' )
for k in range(3,15,3):
ts, xs = midpoint(f0, 0, 1, 1, k)
plt.plot(ts, xs, label=f'{k} steps')
plt.legend()
<matplotlib.legend.Legend at 0x211243d2e70>
To make the picture clearer we can plot the error $x-e^t$ instead of the actual $x$ values.
plt.figure(figsize=(20, 6))
ts = np.linspace(0, 1, 101)
for k in range(3,24,3):
ts, xs = midpoint(f0, 0, 1, 1, k)
plt.plot(ts, xs - np.exp(ts), label=f'{k} steps')
plt.legend()
<matplotlib.legend.Legend at 0x21121dc65d0>
We next give an example where the limitations of Euler's method become very apparent. We want to solve the equation $d^2x/dt^2=-x$ for simple harmonic motion, with $x=1$ and $dx/dt=0$ when $t=0$. As is well known, the exact solution is $x=\cos(t)$ (with velocity $v=dx/dt=-\sin(t)$). In particular, both $x$ and $v$ lie in the interval $[0,1]$.
We cannot use the Euler method or the midpoint method directly for this problem, because it involves second-order derivatives. However, there is a standard trick to deal with this. First note that $dv/dt=d^2x/dt^2=-x$. Consider the vector $u=\begin{bmatrix} x\\v\end{bmatrix}$, and note that
$$ \frac{du}{dt} = \begin{bmatrix} dx/dt\\dv/dt\end{bmatrix} = \begin{bmatrix} v\\ -x\end{bmatrix} = \begin{bmatrix} u[1]\\ -u[0]\end{bmatrix}. $$
This is a first order equation for the vector-valued function $u$, which can be solved by our previous methods. The required values of $x$ are then u[0,0], u[1,0], u[2,0], ..., and we can express this as u[:,0] or (u.T)[0]
We find that the midpoint method works well and gives a result that stays between $-1$ and $1$ as expected, but the errors in the Euler method accumulate making the solution oscillate wildly.
g = lambda t, u: np.array([u[1], -u[0]])
te, ue = euler(g, 0, 20, np.array([1, 0]), 60)
tm, um = midpoint(g, 0, 20, np.array([1, 0]), 60)
xe = ue[:,0] # same as ue.T[0]
xm = (um.T)[0] # same as um[:,0]
plt.figure(figsize=(20, 6))
plt.plot(te, xe, color = 'red', label = 'euler' )
plt.plot(tm, xm, color = 'blue', label = 'midpoint' )
plt.legend()
<matplotlib.legend.Legend at 0x211222795e0>
