Software 3.6.1 (Detecting blow-up). A terminal event stops the integration when \(x\) reaches \(10^8\), which approximates the blow-up time. The bounds of Example 3.4.2 follow:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> def blowup(f, x0, big=1e8, tmax=10):   # first time with x = big
...     ev = lambda t, x: x[0] - big; ev.terminal = True
...     s = solve_ivp(f, (0, tmax), [x0], events=ev, rtol=1e-12, atol=1e-12)
...     return s.t_events[0][0]
...
>>> blowup(lambda t, x: x**2, 1.0), blowup(lambda t, x: x**2, 0.25)
(0.999999989999633, 3.999999989997211)
>>> T = blowup(lambda t, x: t**2 + x**2, 0.0); T
2.003147349427506
>>> s1 = solve_ivp(lambda t, x: t**2 + x**2, (0, 1), [0.0],
...                rtol=1e-12, atol=1e-14)
>>> x1 = s1.y[0, -1]
>>> np.sqrt(np.pi / 2), T, 1 + np.pi / 2 - np.arctan(x1)    # lower, T, upper
(1.2533141373155001, 2.003147349427506, 2.2339149795198647)

Reading the output. For \(x' = x^2\) the computed times are \(1\) and \(4\), matching \(1/x_0\) (the tiny shortfall is the time from \(x = 10^8\) to infinity, \(10^{-8}\)). For \(x' = t^2 + x^2\) the blow-up time is \(T = 2.0031\), inside the comparison bounds \(1.2533\) and \(2.2339\). The integration needs tight tolerances, because the solution grows so fast near \(T\) that the step size must shrink to zero.

MATLAB/Octave: opts = odeset(’Events’, @(t, x) deal(x - 1e8, 1, 0)); then [t, x, te] = ode45(@(t, x) t^2 + x^2, [0 3], 0, opts).

Software 3.6.2 (The variational equation). For the pendulum \(\theta '' = -\sin \theta \) from \(\theta (0) = 1\), \(\theta '(0) = 0\), integrate the system and its variational equation \(Y' = J(t)Y\), with \(J\) the Jacobian of the right-hand side, and compare the first column of \(Y(10)\) with a finite-difference derivative of the flow:

>>> def pend(t, z):                        # pendulum and variational system
...     x, v = z[:2]; Y = z[2:].reshape(2, 2)
...     J = np.array([[0.0, 1.0], [-np.cos(x), 0.0]])
...     return np.concatenate([[v, -np.sin(x)], (J @ Y).ravel()])
...
>>> z0 = np.concatenate([[1.0, 0.0], np.eye(2).ravel()])
>>> Y = solve_ivp(pend, (0, 10), z0, rtol=1e-11, atol=1e-12).y[2:, -1]
>>> Y = Y.reshape(2, 2); np.round(Y, 5)
array([[-0.94351,  0.04995],
       [-1.14356, -0.99933]])
>>> flow = lambda x0: solve_ivp(lambda t, z: [z[1], -np.sin(z[0])], (0, 10),
...                             x0, rtol=1e-11, atol=1e-12).y[:, -1]
>>> h = 1e-6; fd = (flow([1.0 + h, 0.0]) - flow([1.0 - h, 0.0])) / (2 * h)
>>> np.round(fd, 5), round(np.linalg.det(Y), 8)
(array([-0.94351, -1.14356]), 1.0)

Reading the output. The first column of \(Y(10)\), \((-0.94351, -1.14356)\), agrees to five decimals with the finite-difference derivative with respect to \(\theta (0)\), as Theorem 3.5.2 says it must. The determinant of \(Y(10)\) is exactly \(1\): the trace of \(J\) is zero, and Liouville’s formula of the next section then gives \(\det Y = 1\), the statement that the pendulum’s flow preserves area in the phase plane. One integration of the enlarged system gives all the derivatives at once, which is cheaper and more accurate than finite differences.

MATLAB/Octave: append reshape(J*reshape(z(3:6), 2, 2), 4, 1) to the right-hand side and integrate with ode45.

Software 3.6.3 (The flow property). For the logistic equation \(x' = x(1 - x)\), three time units of the flow, and one unit after two:

>>> logi = lambda t, x: x * (1 - x)
>>> phi = lambda T, x0: solve_ivp(logi, (0, T), [x0], rtol=1e-12,
...                               atol=1e-14).y[0, -1]
>>> phi(3.0, 0.1), phi(1.0, phi(2.0, 0.1))  # flow property
(0.6905678577030561, 0.6905678577030469)

Reading the output. \(\phi _3(0.1) = \phi _1(\phi _2(0.1)) = 0.690568\), to \(13\) decimals, illustrating Theorem 3.5.1. The exact value is \(1/(1 + 9e^{-3}) = 0.690568\), from the formula of Example 1.3.3.

MATLAB/Octave: deval(ode45(@(t, x) x.*(1 - x), [0 3], 0.1), 3).