Software 11.6.1 (Energy as a Lyapunov function). The damped pendulum of Example 11.2.4 with \(b = 0.3\), started at rest from \(\theta = 2.5\), close to the top:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> b = 0.3                                 # theta'' + b theta' + sin = 0
>>> pend = lambda t, z: [z[1], -b * z[1] - np.sin(z[0])]
>>> E = lambda z: 0.5 * z[1]**2 + 1 - np.cos(z[0])     # energy
>>> s = solve_ivp(pend, (0, 40), [2.5, 0.0], rtol=1e-10, atol=1e-12,
...               dense_output=True)
>>> ts = np.array([0, 5, 10, 20, 40])
>>> np.round([E(s.sol(t)) for t in ts], 6)
array([1.801144e+00, 6.036960e-01, 1.122340e-01, 5.355000e-03,
       1.600000e-05])
>>> dE = np.diff([E(s.sol(t)) for t in np.linspace(0, 40, 4001)])
>>> dE.max() <= 1e-12                      # E never increases
True
>>> np.round(s.sol(40), 6)
array([ 0.005575, -0.001242])

Reading the output. The energy falls from \(1.80\) to \(1.6 \times 10^{-5}\) in \(40\) time units and never increases at any of \(4000\) sample times, as \(\dot E = -bv^2 \leq 0\) requires. At \(t = 40\) the state is within \(0.006\) of the origin. The decay is not monotone in the angle or the velocity separately, only in the energy: LaSalle’s theorem uses the energy because it is the quantity that friction controls.

MATLAB/Octave: evaluate 0.5*v.^2 + 1 - cos(theta) along the ode45 solution.

Software 11.6.2 (Solving the Lyapunov equation). For \(A = \begin {bmatrix} 0 & 1\\ -2 & -0.5\end {bmatrix}\), solve \(A^{\mathsf T}P + PA = -I\):

>>> from scipy.linalg import solve_continuous_lyapunov
>>> A = np.array([[0.0, 1.0], [-2.0, -0.5]])
>>> np.linalg.eigvals(A)                    # stable spiral
array([-0.25+1.39194109j, -0.25-1.39194109j])
>>> P = solve_continuous_lyapunov(A.T, -np.eye(2))     # A'P + PA = -I
>>> np.round(P, 4), np.round(np.linalg.eigvalsh(P), 4)
(array([[3.125, 0.25 ],
       [0.25 , 1.5  ]]), array([1.4624, 3.1626]))
>>> np.round(A.T @ P + P @ A, 10)
array([[-1.,  0.],
       [ 0., -1.]])

Reading the output. The eigenvalues \(-0.25 \pm 1.39i\) have negative real parts, so by Theorem 11.3.1 the solution \(P\) is positive definite: its eigenvalues are \(1.46\) and \(3.16\). Then \(V = x^{\mathsf T}Px\) satisfies \(\dot V = -\|x\|^2\), which the last line confirms. Note the argument order: SciPy’s routine solves \(MX + XM^{\mathsf H} = C\), so A.T is passed to obtain \(A^{\mathsf T}P + PA\).

MATLAB/Octave: P = lyap(A’, eye(2)) with the Control System Toolbox.

Software 11.6.3 (Explicit and symplectic Euler). The undamped pendulum \(q' = p\), \(p' = -\sin q\) from \((q, p) = (1, 0)\), with step \(h = 0.1\), by explicit Euler, which updates \(q\) and \(p\) from the old values, and symplectic Euler, which updates \(p\) first and uses the new \(p\) to update \(q\):

>>> def euler(h, n, z, symplectic):
...     q, p = z
...     for _ in range(n):
...         if symplectic: p = p - h * np.sin(q); q = q + h * p
...         else: q, p = q + h * p, p - h * np.sin(q)
...     return 0.5 * p**2 + 1 - np.cos(q)
...
>>> 1 - np.cos(1.0)                        # the exact energy H0
0.45969769413186023
>>> for n in (100, 1000, 10000):           # step h = 0.1, from (q, p) = (1, 0)
...     print(n, round(euler(0.1, n, (1.0, 0.0), False), 4),
...           round(euler(0.1, n, (1.0, 0.0), True), 4))
...
100 0.9994 0.4614
1000 5.5392 0.4754
10000 49.0199 0.4727

Reading the output. The exact energy is \(1 - \cos 1 = 0.4597\). Explicit Euler multiplies the area of each small region of phase space by \(1 + h^2\cos q\) at every step, and the energy grows without bound: \(1.0\) after \(100\) steps, \(5.5\) after \(1000\), \(49\) after \(10\,000\), by which time the pendulum is spinning. Symplectic Euler, with the same cost and the same first-order accuracy, keeps the energy within \(3\%\) of its true value however long it runs, because it preserves area exactly. For long integrations of conservative systems, such as planetary orbits and molecular dynamics, the structure of the method matters more than its order.

MATLAB/Octave: the same loops, for k = 1:n, p = p - h*sin(q); q = q + h*p; end.