Software 14.7.1 (The Duffing frequency). The frequency of \(x'' + x + \varepsilon x^3 = 0\) from \(x(0) = 1\), measured from the time between successive maxima, against the Poincaré–Lindstedt value \(1 + \frac 38\varepsilon \):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> def duffing_freq(eps, a=1.0):          # x'' + x + eps x^3 = 0, x(0) = a
...     f = lambda t, z: [z[1], -z[0] - eps * z[0]**3]
...     up = lambda t, z: z[1]; up.direction = -1     # successive maxima
...     s = solve_ivp(f, (0, 60), [a, 0], events=up, rtol=1e-12, atol=1e-12)
...     T = np.diff(s.t_events[0]).mean()
...     return 2 * np.pi / T
...
>>> for eps in (0.01, 0.1, 0.3):
...     print(eps, round(duffing_freq(eps), 6), round(1 + 3 * eps / 8, 6))
...
0.01 1.003742 1.00375
0.1 1.036717 1.0375
0.3 1.106018 1.1125

Reading the output. At \(\varepsilon = 0.01\) the formula is accurate to \(8 \times 10^{-6}\), and at \(\varepsilon = 0.1\) to \(8 \times 10^{-4}\): the error is of order \(\varepsilon ^2\), as the next term of the expansion, \(-\frac {21}{256}\varepsilon ^2a^4\), predicts (\(-0.00082\) at \(\varepsilon = 0.1\)). At \(\varepsilon = 0.3\) the first-order formula is off by \(0.0065\), and higher terms are needed.

MATLAB/Octave: ode45 with an Events function for the maxima.

Software 14.7.2 (Orders of convergence). Euler and fourth-order Runge–Kutta for \(y' = y\), \(y(0) = 1\) on \([0, 1]\), with the step halved three times; the printed numbers are \(\log _2\) of successive error ratios, which estimate the order:

>>> def euler(h):                          # y' = y, y(0) = 1, error at t = 1
...     y = 1.0
...     for _ in range(int(round(1 / h))): y = y + h * y
...     return abs(y - np.e)
...
>>> def rk4(h):
...     y, f = 1.0, (lambda y: y)
...     for _ in range(int(round(1 / h))):
...         k1 = f(y); k2 = f(y + h * k1 / 2); k3 = f(y + h * k2 / 2)
...         k4 = f(y + h * k3); y = y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
...     return abs(y - np.e)
...
>>> hs = [0.1, 0.05, 0.025, 0.0125]
>>> for m in (euler, rk4):
...     e = [m(h) for h in hs]
...     print(m.__name__, np.round(np.log2(np.array(e[:-1]) / e[1:]), 3))
...
euler [0.938 0.968 0.984]
rk4 [3.94  3.97  3.985]

Reading the output. Halving the step halves Euler’s error (observed order \(0.94, 0.97, 0.98 \to 1\)) and divides RK4’s by \(16\) (order \(3.94, 3.97, 3.99 \to 4\)), as Theorem 14.5.1 predicts. For the same accuracy RK4 needs far fewer steps, which more than repays its four evaluations per step.

MATLAB/Octave: the same loops in a script.

Software 14.7.3 (A stiff problem). The equation \(y' = -1000(y - \cos t)\) relaxes in about \(0.001\) time units to its slow solution \(y_p = (10^6\cos t + 10^3\sin t)/(10^6 + 1)\), which varies on a time scale of \(1\). Solved over \([0, 10]\) by an explicit and two implicit methods:

>>> stiff = lambda t, y: [-1000 * (y[0] - np.cos(t))]
>>> yp = lambda t: (1e6 * np.cos(t) + 1e3 * np.sin(t)) / (1e6 + 1)  # slow
>>> for method in ('RK45', 'Radau', 'BDF'):
...     s = solve_ivp(stiff, (0, 10), [0.0], method=method, rtol=1e-6,
...                   atol=1e-9)
...     print(method, s.nfev, abs(s.y[0, -1] - yp(10)))
...
RK45 21548 4.547505373375671e-07
Radau 624 2.0914094900348346e-08
BDF 478 2.0696290237154358e-09

Reading the output. The explicit RK45 needs \(21\,548\) evaluations of \(f\), because stability limits its step to about \(0.003\), and its error at \(t = 10\) is \(4.5 \times 10^{-7}\). The implicit Radau and BDF methods need \(624\) and \(478\) evaluations, a thirty- to forty-fivefold saving, and are more accurate. For non-stiff problems the ranking reverses, and explicit methods are cheaper.

MATLAB/Octave: ode15s (or ode23s) for stiff problems, ode45 otherwise.