Software 12.6.1 (Van der Pol’s cycle from small to large \(\mu \)). A Poincaré section: the event \(x' = 0\) crossed downwards marks each maximum of \(x\). After transients die out, the last maximum is the amplitude of the cycle and the time between the last two maxima its period. The last column is the relaxation approximation \((3 - 2\log 2)\mu \):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> def vdp_cycle(mu, T):                  # x'' - mu (1 - x^2) x' + x = 0
...     f = lambda t, z: [z[1], mu * (1 - z[0]**2) * z[1] - z[0]]
...     top = lambda t, z: z[1]; top.direction = -1   # v = 0 downward: x max
...     s = solve_ivp(f, (0, T), [0.5, 0], events=top, method='LSODA',
...                   rtol=1e-10, atol=1e-12)
...     te, ze = s.t_events[0], s.y_events[0]
...     return round(ze[-1, 0], 4), round(te[-1] - te[-2], 4)
...
>>> for mu, T in [(0.1, 600), (1, 200), (5, 200), (10, 300)]:
...     print(mu, vdp_cycle(mu, T), round((3 - 2 * np.log(2)) * mu, 3))
...
0.1 (2.0001, 6.2871) 0.161
1 (2.0086, 6.6633) 1.614
5 (2.0215, 11.6122) 8.069
10 (2.0143, 19.0784) 16.137

Reading the output. The amplitude stays close to \(2\) for every \(\mu \), as the averaging argument predicts for small \(\mu \) and the relaxation picture for large \(\mu \). The period rises from \(2\pi \) at \(\mu = 0.1\) (\(6.2871\)) through \(6.66\) at \(\mu = 1\) to \(19.08\) at \(\mu = 10\). The leading relaxation formula gives \(16.14\) at \(\mu = 10\); the next term of the asymptotic expansion, \(7.014\mu ^{-1/3} = 3.26\), brings it to \(19.40\), close to the computed value. The stiff solver LSODA is used because for large \(\mu \) the fast jumps and slow drifts have very different time scales.

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

Software 12.6.2 (The Floquet multiplier of the cycle). For \(\mu = 1\), find a point on the cycle and its period from the Poincaré section, then integrate \(\operatorname {div}f = \mu (1 - x^2)\) around one period:

>>> mu = 1.0
>>> f = lambda t, z: [z[1], mu * (1 - z[0]**2) * z[1] - z[0]]
>>> top = lambda t, z: z[1]; top.direction = -1
>>> s = solve_ivp(f, (0, 200), [0.5, 0], events=top, rtol=1e-11, atol=1e-12)
>>> z0, T = s.y_events[0][-1], s.t_events[0][-1] - s.t_events[0][-2]
>>> div = lambda z: mu * (1 - z[0]**2)    # divergence of the vector field
>>> g = lambda t, w: [*f(t, w[:2]), div(w[:2])]
>>> I = solve_ivp(g, (0, T), [*z0, 0.0], rtol=1e-11, atol=1e-12).y[2, -1]
>>> T, I, np.exp(I)                        # period, integral, multiplier
(6.663286859314979, -7.058932808723691, 0.0008596950636685128)

Reading the output. The divergence integral over one period is \(-7.06\), so the non-trivial multiplier is \(e^{-7.06} = 0.00086\): a displacement off the cycle is reduced more than a thousand-fold in one revolution. The cycle is strongly attracting, which is why every orbit in Figure 16 joins it within a few turns. The integral is negative although \(\operatorname {div}f > 0\) for \(|x| < 1\): the orbit spends enough of its time at \(|x| > 1\).

MATLAB/Octave: integrate the augmented system [f; mu*(1 - x^2)] over one period.

Software 12.6.3 (A cycle with no formula). For Example 12.2.3, two solutions, one starting near the origin and one far outside, and their radius at late times:

>>> h = lambda t, z: [z[0] - z[1] - z[0]**3, z[0] + z[1] - z[1]**3]
>>> r = lambda z: np.hypot(*z)
>>> for z0 in ([0.01, 0.0], [3.0, 3.0]):
...     s = solve_ivp(h, (0, 60), z0, rtol=1e-10, atol=1e-12,
...                   dense_output=True)
...     print(z0, [round(r(s.sol(t)), 4) for t in (40, 50, 60)])
...
[0.01, 0.0] [1.2478, 1.193, 1.0972]
[3.0, 3.0] [1.2461, 1.1821, 1.0894]

Reading the output. Both radii lie between \(1\) and \(\sqrt 2 = 1.414\) at late times, as the trapping annulus requires, and at each time the two values almost coincide: the solutions have converged to the same periodic orbit, along which the radius oscillates between about \(1.09\) and \(1.25\) because the cycle is not a circle. The equilibrium check of Example 12.2.3 (only the origin) was done by computing a resultant of the two right-hand sides, whose only real root is \(x = 0\).

MATLAB/Octave: ode45 and hypot on the solution.