Software 7.6.1 (Monodromy matrices and multipliers). The monodromy matrix of Mathieu’s equation is computed by integrating two solutions over one period; its columns are the solutions with initial values \((1, 0)\) and \((0, 1)\):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> from scipy.optimize import brentq
>>> def monodromy(d, e, T=2 * np.pi):     # y'' + (d + e cos t) y = 0
...     f = lambda t, z: [z[1], -(d + e * np.cos(t)) * z[0],
...                       z[3], -(d + e * np.cos(t)) * z[2]]
...     z = solve_ivp(f, (0, T), [1, 0, 0, 1], rtol=1e-11,
...                   atol=1e-12).y[:, -1]
...     return np.array([[z[0], z[2]], [z[1], z[3]]])   # columns: solutions
...
>>> for d, e in [(0.5, 0.2), (0.25, 0.2), (1.0, 0.2)]:
...     M = monodromy(d, e); rho = np.linalg.eigvals(M)
...     print(d, e, round(np.trace(M), 4), np.round(np.abs(rho), 4),
...           round(np.linalg.det(M), 8))
...
0.5 0.2 -0.7026 [1. 1.] 1.0
0.25 0.2 -2.393 [0.5395 1.8534] 1.0
1.0 0.2 2.0005 [0.9769 1.0237] 1.0

Reading the output. At \((\delta , \varepsilon ) = (0.5, 0.2)\), between the first two tongues, \(|\Delta | = 0.70 < 2\) and both multipliers lie on the unit circle: stable. At \((0.25, 0.2)\), in the middle of the first tongue, \(\Delta = -2.39\) and the multipliers are real, \(-1.853\) and \(-0.540\) (moduli shown), so a solution grows by the factor \(1.853\), changing sign, every period: parametric resonance at half the forcing frequency. At \((1, 0.2)\), \(\Delta = 2.0005\) lies just inside the narrow second tongue. The determinant is \(1\) in every case, as Theorem 7.3.2(i) requires for Hill’s equation.

MATLAB/Octave: integrate the four-dimensional system with ode45 over \([0, 2\pi ]\) and eig the resulting \(2 \times 2\) matrix.

Software 7.6.2 (The edges of the first tongue). The boundaries of the first tongue at \(\varepsilon = 0.2\) are the roots of \(\Delta (\delta ) = -2\), found by brentq, and compared with harmonic balance:

>>> tr = lambda d, e: np.trace(monodromy(d, e))
>>> e = 0.2
>>> lo = brentq(lambda d: tr(d, e) + 2, 0.10, 0.25)   # first tongue edges
>>> hi = brentq(lambda d: tr(d, e) + 2, 0.25, 0.40)
>>> (lo, hi), (0.25 - e / 2, 0.25 + e / 2)
((0.14524515179353695, 0.34474668423949), (0.15, 0.35))

Reading the output. The numerical edges, \(0.1452\) and \(0.3447\), are close to the first-order values \(0.15\) and \(0.35\). The small shift of both edges to the left is the second-order correction: the expansion of the tongue boundaries continues as \(\delta = \frac 14 \pm \frac \varepsilon 2 - \frac {\varepsilon ^2}{8} + O(\varepsilon ^3)\), and \(-\varepsilon ^2/8 = -0.005\) gives \(0.145\) and \(0.345\). Because the tongue boundary is where a \(4\pi \)-periodic solution exists, \(\Delta + 2\) changes sign there, which makes root finding reliable.

MATLAB/Octave: fzero(@(d) trace(monodromy(d, 0.2)) + 2, [0.1 0.25]).

Software 7.6.3 (Kapitza’s pendulum). The full non-linear equation \(\theta '' + (\delta + \varepsilon \cos \tau )\sin \theta = 0\), with \(\theta \) measured from the upright position, started at \(\theta = 0.1\) and integrated for \(200\) periods, on both sides of the boundary \(\delta = -\varepsilon ^2/2 = -0.045\) for \(\varepsilon = 0.3\):

>>> def pendulum(d, e, theta0=0.1, periods=200):    # vibrated support
...     f = lambda t, z: [z[1], -(d + e * np.cos(t)) * np.sin(z[0])]
...     s = solve_ivp(f, (0, 2 * np.pi * periods), [theta0, 0.0],
...                   rtol=1e-9, atol=1e-11, max_step=0.05)
...     return np.abs(s.y[0]).max()
...
>>> [(d, round(pendulum(d, 0.3), 3)) for d in (-0.03, -0.06)]
[(-0.03, 0.1), (-0.06, 966.889)]

Reading the output. With \(\delta = -0.03\), inside the stable strip, the largest deviation in \(200\) periods is the initial \(0.1\): the vibration holds the pendulum upright. With \(\delta = -0.06\), outside it, the pendulum falls over and goes on to rotate, the angle reaching \(967\) radians. The linear Floquet analysis predicts the non-linear behaviour near the equilibrium, as the stability theory of the following sections will justify. In physical terms, for a pendulum of length \(\ell = 0.2\) m shaken with amplitude \(a = 1\) cm, stability needs \(\omega ^2 > 2g\ell /a^2\), that is \(\omega > 198\) rad/s, about \(32\) shakes a second (Problem 7.7.5).

MATLAB/Octave: ode45(@(t, z) [z(2); -(d + e*cos(t))*sin(z(1))], [0 400*pi], [0.1; 0]).