Software 13.5.1 (A Hopf bifurcation in the Brusselator). With \(a = 1\), so \(b_c = 2\), the half-height of the oscillation of \(x\) after transients, and its ratio to \(\sqrt {b - b_c}\):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> a = 1.0                                 # Brusselator: Hopf at b = 2
>>> def amplitude(b):
...     f = lambda t, z: [a - (b + 1) * z[0] + z[0]**2 * z[1],
...                       b * z[0] - z[0]**2 * z[1]]
...     s = solve_ivp(f, (0, 600), [a + 0.1, b / a], rtol=1e-9, atol=1e-11,
...                   dense_output=True)
...     x = s.sol(np.linspace(500, 600, 20001))[0]
...     return (x.max() - x.min()) / 2
...
>>> for b in (1.9, 2.02, 2.08, 2.18, 2.32):
...     A = amplitude(b)
...     print(b, round(A, 4), round(A / np.sqrt(max(b - 2, 1e-12)), 3))
...
1.9 0.0 0.013
2.02 0.1648 1.165
2.08 0.3385 1.197
2.18 0.5306 1.251
2.32 0.7511 1.328

Reading the output. Below the bifurcation (\(b = 1.9\)) the oscillation dies out and the amplitude is \(0\). Above it, the amplitude grows from \(0.16\) at \(b = 2.02\) to \(0.75\) at \(b = 2.32\), and its ratio to \(\sqrt {b - 2}\) is nearly constant close to the bifurcation, about \(1.17\), rising slowly further away as higher-order terms enter: the square-root law of Theorem 13.3.1(i) for a supercritical Hopf bifurcation.

MATLAB/Octave: sweep \(b\) with ode45 and take max and min of \(x\) over the last part of each solution.

Software 13.5.2 (Sensitive dependence and the Lyapunov exponent). Two orbits of the Lorenz system starting \(10^{-9}\) apart on the attractor, and the largest Lyapunov exponent estimated by repeatedly renormalising the separation of two nearby orbits (Benettin’s method):

>>> sig, rho, beta = 10.0, 28.0, 8 / 3
>>> lor = lambda t, u: [sig * (u[1] - u[0]), u[0] * (rho - u[2]) - u[1],
...                     u[0] * u[1] - beta * u[2]]
>>> u0 = solve_ivp(lor, (0, 50), [1, 1, 1], rtol=1e-10).y[:, -1]  # attractor
>>> ts = np.linspace(0, 25, 6)
>>> A = solve_ivp(lor, (0, 25), u0, t_eval=ts, rtol=1e-12, atol=1e-12).y
>>> B = solve_ivp(lor, (0, 25), u0 + [1e-9, 0, 0], t_eval=ts, rtol=1e-12,
...               atol=1e-12).y
>>> np.round(np.linalg.norm(A - B, axis=0), 9)
array([1.00000000e-09, 1.05000000e-07, 4.81600000e-06, 2.11337590e-02,
       1.34042253e-01, 3.13333813e+00])
>>> d0, total, u = 1e-8, 0.0, u0.copy(); v = u + [d0, 0, 0]
>>> for k in range(400):                    # renormalise every 0.5 time units
...     u = solve_ivp(lor, (0, 0.5), u, rtol=1e-10, atol=1e-12).y[:, -1]
...     v = solve_ivp(lor, (0, 0.5), v, rtol=1e-10, atol=1e-12).y[:, -1]
...     d = np.linalg.norm(v - u); total += np.log(d / d0)
...     v = u + (v - u) * d0 / d
...
>>> total / (400 * 0.5)                     # largest Lyapunov exponent
0.9003744690325769

Reading the output. The separation grows from \(10^{-9}\) to \(10^{-7}\) in \(5\) time units, to \(5 \times 10^{-6}\) by \(10\), and to \(3\), the size of the attractor, by \(25\). The exponent estimated over \(200\) time units is \(0.900\), against the accepted value \(0.9056\); the estimate converges slowly, like the square root of the averaging time. The separation does not grow exactly exponentially over short intervals, because the stretching rate varies around the attractor, and the exponent is its long-time average.

MATLAB/Octave: the same renormalisation loop with ode45.

Software 13.5.3 (Period doubling in the logistic map). The number of distinct points visited after transients, for five values of \(r\), and the Feigenbaum ratios from the first five doubling values:

>>> def attractor(r, n=2000, keep=64):
...     x = 0.5
...     for _ in range(n): x = r * x * (1 - x)
...     pts = set()
...     for _ in range(keep):
...         x = r * x * (1 - x); pts.add(round(x, 6))
...     return sorted(pts)
...
>>> [(r, len(attractor(r))) for r in (2.8, 3.2, 3.5, 3.56, 3.9)]
[(2.8, 1), (3.2, 2), (3.5, 4), (3.56, 8), (3.9, 64)]
>>> r = [3.0, 1 + np.sqrt(6), 3.5440903, 3.5644073, 3.5687594]  # doublings
>>> [round((r[k] - r[k-1]) / (r[k+1] - r[k]), 3) for k in (1, 2, 3)]
[4.751, 4.656, 4.668]

Reading the output. The attractor has \(1\), \(2\), \(4\) and \(8\) points at \(r = 2.8\), \(3.2\), \(3.5\) and \(3.56\), and at \(r = 3.9\) every one of the \(64\) sampled points is different: chaos. The ratios of successive intervals between doublings, \(4.751\), \(4.656\), \(4.668\), approach Feigenbaum’s \(\delta = 4.6692\).

MATLAB/Octave: a for loop with unique(round(x, 6)).