Software 4.6.1 (Liouville’s formula). For \(x'' + 0.2x' + (1 + 0.5\cos t)x = 0\), a damped oscillator with a periodically varying stiffness, the fundamental matrix is computed by integrating \(\Phi ' = A(t)\Phi \) from \(\Phi (0) = I\); here \(\tr A = -0.2\):

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> A = lambda t: np.array([[0.0, 1.0], [-(1 + 0.5 * np.cos(t)), -0.2]])
>>> def fund(t, z):                        # Phi' = A(t) Phi, Phi(0) = I
...     return (A(t) @ z.reshape(2, 2)).ravel()
...
>>> sol = solve_ivp(fund, (0, 20), np.eye(2).ravel(), rtol=1e-11, atol=1e-13,
...                 dense_output=True)
>>> for t in (5.0, 10.0, 20.0):
...     Phi = sol.sol(t).reshape(2, 2)
...     print(t, round(np.linalg.det(Phi), 8), round(np.exp(-0.2 * t), 8))
...
5.0 0.36787944 0.36787944
10.0 0.13533528 0.13533528
20.0 0.01831564 0.01831564

Reading the output. The determinant of the numerically computed \(\Phi (t)\) equals \(e^{-0.2t}\) to eight decimals at every time, although \(\Phi \) itself has no closed form: the stiffness varies with time, but only the trace matters for the volume. This is a useful check on numerical work: a computed fundamental matrix whose determinant drifts from Liouville’s value has accumulated error.

MATLAB/Octave: integrate @(t, z) reshape(A(t)*reshape(z, 2, 2), 4, 1) with ode45 from reshape(eye(2), 4, 1).

Software 4.6.2 (Variation of parameters). For \(y'' + y = \sec t\) on \((-\pi /2, \pi /2)\), with \(y_1 = \cos t\), \(y_2 = \sin t\):

>>> import sympy as sp
>>> t = sp.symbols('t')
>>> y1, y2, g = sp.cos(t), sp.sin(t), 1 / sp.cos(t)   # y'' + y = sec t
>>> W = sp.simplify(y1 * y2.diff(t) - y2 * y1.diff(t)); W
1
>>> yp = sp.simplify(-y1 * sp.integrate(y2 * g / W, t)
...                  + y2 * sp.integrate(y1 * g / W, t)); yp
t*sin(t) + log(cos(t))*cos(t)
>>> sp.simplify(yp.diff(t, 2) + yp - g)
0

Reading the output. The Wronskian is \(1\), and the particular solution is \(y_p = t\sin t + \cos t\log \cos t\), which the last line confirms by substitution. The method of undetermined coefficients cannot produce this, because \(\sec t\) is not of the exponential-polynomial form it requires.

MATLAB/Octave: syms y(t); dsolve(diff(y, 2) + y == sec(t)).

Software 4.6.3 (A non-commuting coefficient matrix). For \(B(t) = \begin {bmatrix} 1 & t\\ 0 & 0\end {bmatrix}\), compare the true transition matrix \(\Phi (1, 0)\) with \(\exp \int _0^1B(s)\,ds\):

>>> from scipy.linalg import expm
>>> B = lambda t: np.array([[1.0, t], [0.0, 0.0]])
>>> intB = lambda t: np.array([[t, t**2 / 2], [0.0, 0.0]])
>>> B(1.0) @ intB(1.0) - intB(1.0) @ B(1.0)        # they do not commute
array([[ 0. , -0.5],
       [ 0. ,  0. ]])
>>> Phi = solve_ivp(lambda t, z: (B(t) @ z.reshape(2, 2)).ravel(), (0, 1),
...                 np.eye(2).ravel(), rtol=1e-12, atol=1e-13).y[:, -1]
>>> np.round(Phi.reshape(2, 2), 6), np.round(expm(intB(1.0)), 6)
(array([[2.718282, 0.718282],
       [0.      , 1.      ]]), array([[2.718282, 0.859141],
       [0.      , 1.      ]]))

Reading the output. \(B(1)\) and \(\int _0^1B\) do not commute, and the two matrices differ in their top-right entry: the true value is \(e - 2 = 0.718282\), obtained by solving \(x_2 = c\), \(x_1' = x_1 + tc\) directly, while the exponential gives \((e - 1)/2 = 0.859141\). Proposition 4.5.1 cannot be used without checking commutativity.

MATLAB/Octave: expm([1 0.5; 0 0]) for the exponential.