Software 5.6.1 (Jordan form and the exponential in SymPy). The defective matrix of Example 5.2.2, and a \(3 \times 3\) matrix with a complex pair:

>>> import sympy as sp
>>> t = sp.symbols('t', real=True)
>>> A = sp.Matrix([[1, 1], [-1, 3]])
>>> A.eigenvals()                          # one eigenvalue, multiplicity 2
{2: 2}
>>> P, J = A.jordan_form(); J
Matrix([
[2, 1],
[0, 2]])
>>> sp.simplify((A * t).exp())
Matrix([
[(1 - t)*exp(2*t),       t*exp(2*t)],
[     -t*exp(2*t), (t + 1)*exp(2*t)]])
>>> B = sp.Matrix([[0, -2, 0], [1, -2, 0], [0, 0, -1]])
>>> B.eigenvals()
{-1 - I: 1, -1 + I: 1, -1: 1}
>>> sp.simplify((B * t).exp())[0, :]
Matrix([[sqrt(2)*exp(-t)*sin(t + pi/4), -2*exp(-t)*sin(t), 0]])

Reading the output. SymPy finds the double eigenvalue \(2\) and the single Jordan block, and \(e^{At}\) agrees with Example 5.2.2; the factor \(t\) in every entry is the signature of the missing eigenvector. For \(B\), the eigenvalues are \(-1\) and \(-1 \pm i\), and the first row of \(e^{Bt}\) contains \(e^{-t}\sin (t + \pi /4)\) and \(e^{-t}\sin t\), decaying oscillations at the frequency given by the imaginary part, as Theorem 5.2.1(iii) predicts.

MATLAB/Octave: [V, J] = jordan(sym([1 1; -1 3])) (Symbolic Math Toolbox) and expm(A*t) numerically.

Software 5.6.2 (Classifying planar systems). A direct implementation of Theorem 5.4.1, applied to five matrices:

>>> import numpy as np
>>> def classify(A):
...     tau, delta = np.trace(A), np.linalg.det(A)
...     disc = tau**2 - 4 * delta
...     s = 'stable' if tau < 0 else 'unstable'
...     if delta < 0: kind = 'saddle'
...     elif delta == 0: kind = 'degenerate (zero eigenvalue)'
...     elif tau == 0: kind = 'centre'
...     else: kind = s + (' spiral' if disc < 0 else ' node')
...     return round(tau, 3), round(delta, 3), kind
...
>>> for M in ([[0, 1], [-4, -1]], [[2, 1], [3, -1]], [[-3, 1], [1, -3]],
...           [[0, 2], [-2, 0]], [[1, -5], [5, 1]]):
...     print(M, classify(np.array(M, float)))
...
[[0, 1], [-4, -1]] (-1.0, 4.0, 'stable spiral')
[[2, 1], [3, -1]] (1.0, -5.0, 'saddle')
[[-3, 1], [1, -3]] (-6.0, 8.0, 'stable node')
[[0, 2], [-2, 0]] (0.0, 4.0, 'centre')
[[1, -5], [5, 1]] (2.0, 26.0, 'unstable spiral')

Reading the output. The first matrix is the system form of the damped oscillator \(y'' + y' + 4y = 0\), with \(\tau = -1\), \(\Delta = 4\) and \(\tau ^2 < 4\Delta \): a stable spiral, a decaying oscillation. The second has \(\Delta < 0\), a saddle. The symmetric third matrix has eigenvalues \(-2\) and \(-4\), a stable node. The fourth is a centre, and the fifth an unstable spiral with eigenvalues \(1 \pm 5i\). In practice, comparing \(\tau ^2\) with \(4\Delta \) in floating point is delicate near the parabola, where the classification is structurally unstable: an arbitrarily small change in \(A\) moves a degenerate node into either a node or a spiral.

MATLAB/Octave: trace(A), det(A) and eig(A).

Software 5.6.3 (Transient growth). The matrix \(N = \begin {bmatrix} -1 & 10\\ 0 & -2\end {bmatrix}\) has eigenvalues \(-1\) and \(-2\), so every solution of \(x' = Nx\) decays. Compute \(\|e^{Nt}\|\) over time:

>>> from scipy.linalg import expm
>>> N = np.array([[-1.0, 10.0], [0.0, -2.0]])
>>> np.linalg.eigvals(N)
array([-1., -2.])
>>> ts = np.linspace(0, 6, 601)
>>> norms = np.array([np.linalg.norm(expm(N * s), 2) for s in ts])
>>> ts[norms.argmax()], norms.max(), norms[-1]
(0.66, 2.563475247123823, 0.02485001537720811)

Reading the output. Although both eigenvalues are negative, \(\|e^{Nt}\|\) rises to \(2.56\) at \(t = 0.66\) before decaying: some initial states are amplified more than two and a half times before the decay predicted by the eigenvalues takes over. The cause is the large off-diagonal entry, which makes the eigenvectors \((1, 0)\) and \((10, -1)\) nearly parallel. By \(t = 6\) the norm is \(0.025\), roughly \(Ke^{-t}\) with \(K \approx 10\). Eigenvalues describe the long run; the constant \(K\) in Theorem 5.3.2(i) describes the short run.

MATLAB/Octave: arrayfun(@(s) norm(expm(N*s)), linspace(0, 6, 601)).