Software 2.6.1 (Picard iteration in SymPy). The iterates of Example 2.3.2, computed by exact integration, compared with the Taylor series of \(\tan t\):

>>> import sympy as sp
>>> t, s = sp.symbols('t s')
>>> f = lambda t, x: 1 + x**2              # x' = 1 + x^2, x(0) = 0
>>> phi = [sp.Integer(0)]
>>> for k in range(4):                     # the Picard iteration
...     p = phi[-1].subs(t, s)
...     phi.append(sp.expand(sp.integrate(f(s, p), (s, 0, t))))
...
>>> phi[3]
t**7/63 + 2*t**5/15 + t**3/3 + t
>>> sp.series(sp.tan(t), t, 0, 9)
t + t**3/3 + 2*t**5/15 + 17*t**7/315 + O(t**9)
>>> [sp.degree(p, t) for p in phi[1:]]
[1, 3, 7, 15]

Reading the output. \(\phi _3\) agrees with \(\tan t = t + \frac {t^3}{3} + \frac {2t^5}{15} + \frac {17t^7}{315} + \cdots \) up to \(t^5\); its \(t^7\) coefficient \(\frac {1}{63} = \frac {5}{315}\) is still wrong. The degrees \(1, 3, 7, 15\) double (plus one) at each step, because squaring a polynomial doubles its degree, so the iterates contain many wrong high-order terms that later steps correct.

MATLAB/Octave: syms t s; p = sym(0); for k = 1:3, p = int(1 + subs(p, t, s)^2, s, 0, t); end.

Software 2.6.2 (Convergence and the guaranteed interval). The errors \(\phi _k - \tan t\) at three times, and the best interval that Theorem 2.3.1 can promise over the choice of \(b\):

>>> import numpy as np
>>> xs = (0.25, 0.5, 1.0)
>>> for k in range(1, 5):
...     g = sp.lambdify(t, phi[k])
...     print(k, [round(float(g(x) - np.tan(x)), 5) for x in xs])
...
1 [-0.00534, -0.0463, -0.55741]
2 [-0.00013, -0.00464, -0.22407]
3 [-0.0, -0.00035, -0.07487]
4 [-0.0, -2e-05, -0.02045]
>>> b = np.linspace(0.01, 5, 2000); alpha = np.minimum(1.0, b / (1 + b**2))
>>> b[np.argmax(alpha)], alpha.max()
(1.0010105052526264, 0.49999974497761523)

Reading the output. At \(t = 0.25\) and \(t = 0.5\), inside the guaranteed interval, the error falls by an order of magnitude or more at each step; at \(t = 1\), outside it, convergence is slower but still visible (\(-0.557\), \(-0.224\), \(-0.075\), \(-0.020\)). The iterates in fact converge on the whole of \((-\pi /2, \pi /2)\); the theorem is a sufficient condition, not a description. The search over \(b\) confirms that \(\alpha = b/(1 + b^2)\) is maximised at \(b = 1\), with \(\alpha = \frac 12\).

MATLAB/Octave: b = linspace(0.01, 5, 2000); [amax, i] = max(min(1, b./(1 + b.^2))).

Software 2.6.3 (Continuous dependence and its failure). Two solutions of \(x' = 2\sin x\) from initial values \(0.001\) apart, against the bound of Theorem 2.4.2; then the equation \(y' = 3y^{2/3}\) of Example 1.3.2 from \(0\) and from \(10^{-12}\):

>>> from scipy.integrate import solve_ivp
>>> L, T = 2.0, 1.0                        # 2 sin x has Lipschitz L = 2
>>> rhs = lambda t, x: 2 * np.sin(x)
>>> x1 = solve_ivp(rhs, (0, T), [1.0], rtol=1e-11, atol=1e-12)
>>> x2 = solve_ivp(rhs, (0, T), [1.001], rtol=1e-11, atol=1e-12)
>>> gap = abs(x2.y[0, -1] - x1.y[0, -1]); gap, 0.001 * np.exp(L * T)
(0.0005542867025911669, 0.007389056098930651)
>>> nu = lambda t, y: 3 * np.cbrt(y)**2    # y' = 3 y^(2/3): uniqueness fails
>>> [solve_ivp(nu, (0, 1), [y0], rtol=1e-10, atol=1e-14).y[0, -1]
...  for y0 in (0.0, 1e-12)]
[0.0, 1.0003001102501183]

Reading the output. After one unit of time the two solutions are \(0.00055\) apart, well within the Gronwall bound \(0.001e^2 = 0.0074\); the bound allows for the worst case, in which every direction is expanding at the maximal rate \(L\), while near \(x = 1\) the solutions are being drawn towards the stable equilibrium \(\pi \). For the non-Lipschitz equation, continuous dependence fails completely: starting at \(0\) the solver returns \(0\), but starting at \(10^{-12}\) it follows \(y \approx t^3\) to \(y(1) = 1.0003\). An initial difference of \(10^{-12}\) has become a difference of \(1\), because the two starting points lie on different members of the infinite family of solutions through the origin.

MATLAB/Octave: [t, x] = ode45(@(t, x) 2*sin(x), [0 1], 1.001).