Ordinary Differential Equations (7): Stability Theory

A small push hits a system. Does it return to rest, drift away, or break entirely? That single question decides whether bridges survive storms, ecosystems recover from droughts, and economies bounce back from crises. Stability theory answers it — and it does so without ever solving the differential equation. We will learn to read the destiny of a system off the geometry of its phase plane.

What You Will Learn#

• Three precise notions: Lyapunov stable, asymptotically stable, unstable
• Linearization via the Jacobian and the Hartman-Grobman theorem
• Lyapunov’s direct method — proving stability with energy-like functions
• LaSalle’s invariance principle for borderline cases
• Trace-determinant classification of all 2D linear systems
• Four canonical bifurcations: saddle-node, transcritical, pitchfork, Hopf
• Worked applications: pendulum, predator-prey, inverted pendulum control

Prerequisites#

• Chapter 6 : linear systems, eigenvalues, phase portraits
• Multivariable calculus: partial derivatives, Jacobian matrix

A Visual Tour Before the Theory#

Stability is, at heart, a geometric statement about how trajectories move in phase space. Six pictures tell the entire story of 2D linear systems.

Reading the six pictures. Each portrait is geometry made out of two eigenvalues $\lambda_{1,2}$ of the linearization. Read them in pairs.

• Saddle (real eigenvalues of opposite sign, $\lambda_1 > 0 > \lambda_2$ ). The two eigenvectors carve phase space into a cross: along the unstable manifold $E^u$ trajectories shoot outward at rate $e^{\lambda_1 t}$ , along the stable manifold $E^s$ they decay inward at rate $e^{\lambda_2 t}$ . Every other trajectory is a hyperbola threading between the two: it approaches the equilibrium along $E^s$ , swerves at the closest point, and leaves along $E^u$ . The pass between two mountains is the canonical picture — uphill in one direction, downhill in the perpendicular one, and you can balance on the saddle for an instant if you arrive on $E^s$ exactly.
• Stable / unstable node (real eigenvalues of the same sign). All trajectories enter (or leave) the equilibrium tangent to the eigenvector with the slower (smaller-magnitude) eigenvalue, because along the fast direction the contribution dies first. The picture is a “comb” of curves combed into the slow eigendirection.
• Stable / unstable spiral (complex pair $\alpha \pm i\beta$ ). The real part $\alpha$ is a radial decay/growth rate; the imaginary part $\beta$ is an angular frequency. In polar coordinates the linearized flow is $\dot r = \alpha r,\ \dot\theta = \beta$ , so trajectories trace logarithmic spirals — distance from the origin shrinks (or grows) like $e^{\alpha t}$ while the angle rotates at constant rate $\beta$ . The direction of rotation (clockwise vs counterclockwise) is fixed by the off-diagonal entries of $A$ , not by the eigenvalues themselves.
• Center (purely imaginary eigenvalues). $\alpha = 0$ kills the radial motion and only rotation survives; orbits are nested ellipses determined by $A$ . Centers are the only neutrally stable case in 2D linear systems and are structurally fragile — any small perturbation that adds dissipation turns the center into a spiral.
• Degenerate node (repeated real eigenvalue with one eigenvector). Only one eigendirection exists, and a generalized eigenvector contributes a $t\,e^{\lambda t}$ term, which bends every trajectory tangent to that eigenvector as it approaches the equilibrium. It is the pinched limit between a node (two eigendirections) and a spiral (none).

The single rule of thumb: the real parts decide whether you decay or grow; the imaginary parts decide whether you rotate; the eigenvectors set the axes the picture is drawn on.

For nonlinear systems, the same six pictures still appear — but only locally, near each equilibrium. The damped pendulum and the Lotka-Volterra predator-prey model both show this beautifully:

Stability Defined Precisely#

Consider $\mathbf{x}' = \mathbf{f}(\mathbf{x})$ with equilibrium $\mathbf{x}^*$ (so $\mathbf{f}(\mathbf{x}^*) = \mathbf{0}$ ).

Intuition: nearby trajectories stay nearby forever.

Asymptotically stable. Lyapunov stable and $\mathbf{x}(t) \to \mathbf{x}^*$ as $t \to \infty$ . Intuition: nearby trajectories not only stay nearby but eventually return.

Unstable. Not Lyapunov stable.

The basin of attraction is the set of all initial conditions that converge to $\mathbf{x}^*$ . Asymptotic stability is a local property; the basin tells you how local.

Why two definitions? A center (closed orbits) is Lyapunov stable but not asymptotically stable — trajectories stay close but never settle. The Lotka-Volterra model is the classic example.

Linearization: The Jacobian Method#

Hartman-Grobman theorem#

If every eigenvalue of $J$ has nonzero real part (a hyperbolic equilibrium), then the nonlinear system is locally topologically equivalent to its linearization $\mathbf{u}' = J\mathbf{u}$ .

• All $\operatorname{Re}(\lambda) < 0$ : asymptotically stable
• Any $\operatorname{Re}(\lambda) > 0$ : unstable
• Purely imaginary eigenvalues: linearization fails — use Lyapunov methods

Example: damped pendulum#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

gamma, omega0 = 0.3, 1.0
def pendulum(X, t):
    x, y = X
    return [y, -gamma*y - omega0**2*np.sin(x)]

fig, ax = plt.subplots(figsize=(10, 6))
t = np.linspace(0, 50, 2000)
for x0 in np.linspace(-3*np.pi, 3*np.pi, 15):
    for y0 in np.linspace(-3, 3, 5):
        sol = odeint(pendulum, [x0, y0], t)
        ax.plot(sol[:,0], sol[:,1], 'b-', linewidth=0.3, alpha=0.5)
for n in range(-2, 3):
    ax.plot(n*np.pi, 0, 'go' if n%2==0 else 'rx', markersize=10)
ax.set(xlabel='theta', ylabel='omega', title='Damped pendulum phase portrait')
ax.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

Lyapunov’s Direct Method#

The big idea#

Stability without solving the ODE. Construct an energy-like scalar function $V(\mathbf{x})$ and show it decreases along trajectories.

Requirements.

1. $V(\mathbf{x}^*) = 0$ and $V(\mathbf{x}) > 0$ otherwise (positive definite)
2. $\dot V = \nabla V \cdot \mathbf{f}(\mathbf{x}) \leq 0$ along orbits

Stability theorems.

Why it works — the picture#

Trajectories cross level sets of $V$ inward. Since $V$ has a minimum at $\mathbf{x}^*$ , they cannot escape arbitrarily-small level sets, and (with strict descent) they slide all the way to the bottom.

How to find $V$ #

• Physical energy: kinetic + potential for mechanical systems
• Quadratic forms: $V = \mathbf{x}^T P \mathbf{x}$ where $P$ solves the Lyapunov equation $A^T P + PA = -Q$
• Trial: start with $V = x^2 + y^2$ , compute $\dot V$ , adjust coefficients

Example: pendulum energy#

The hanging position is stable. LaSalle’s principle (next) upgrades this to asymptotic.

LaSalle’s Invariance Principle#

Sometimes $\dot V \leq 0$ but vanishes on a whole set, not just at $\mathbf{x}^*$ . Standard Lyapunov only gives stability, not attraction.

Theorem. Let $E = \{\mathbf{x} : \dot V = 0\}$ and $M$ be the largest invariant subset of $E$ . Every bounded trajectory approaches $M$ .

If $M = \{\mathbf{x}^*\}$ , then $\mathbf{x}^*$ is asymptotically stable.

For the damped pendulum, $\dot V = 0$ requires $\omega = 0$ . But on the line $\omega = 0$ the dynamics force $\dot\omega = -\omega_0^2 \sin\theta \neq 0$ unless also $\theta = 0$ . So $M = \{(0,0)\}$ and we get asymptotic stability for free.

Trace-Determinant Classification#

For $\mathbf{x}' = A\mathbf{x}$ in 2D, set $\tau = \operatorname{tr}(A)$ and $\Delta = \det(A)$ . Then $\lambda_{1,2} = \tfrac12(\tau \pm \sqrt{\tau^2 - 4\Delta})$ , and:

A single picture compresses this entire table:

Bifurcations: When the Picture Itself Changes#

Slowly turn a parameter knob $r$ . Equilibria can be born, die, or swap stability. Four “normal forms” capture every codimension-1 bifurcation locally.

Hopf is the mechanism behind every self-sustained oscillation in nature — from heartbeats to the pulsing of variable stars.

Application 1: Lotka-Volterra Predator-Prey#

makes every orbit closed. The system has periodic population cycles (right panel of fig 2).

Application 2: Inverted Pendulum Control#

The inverted equilibrium is a saddle. Linear feedback $u = -K\mathbf{x}$ shifts the closed-loop eigenvalues into the open left half-plane, converting the saddle into a stable focus.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import numpy as np
from scipy.linalg import solve_continuous_are
from scipy.integrate import odeint
import matplotlib.pyplot as plt

g, l, m, M_cart = 9.81, 0.5, 0.1, 1.0
A = np.array([[0,1,0,0],[0,0,-m*g/M_cart,0],
              [0,0,0,1],[0,0,(M_cart+m)*g/(M_cart*l),0]])
B = np.array([[0],[1/M_cart],[0],[-1/(M_cart*l)]])

Q = np.diag([10, 1, 100, 1])
R = np.array([[0.1]])
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print("Closed-loop eigenvalues:", np.linalg.eigvals(A - B @ K))

t = np.linspace(0, 5, 500)
sol = odeint(lambda X, t: (A @ X + B @ (-K @ X)).flatten(),
             [0,0,0.2,0], t)

plt.figure(figsize=(8, 4))
plt.plot(t, sol[:,2]*180/np.pi, 'b-', linewidth=2)
plt.xlabel('Time (s)'); plt.ylabel('Angle (degrees)')
plt.title('Inverted pendulum stabilized by LQR')
plt.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

Exercises#

Basic.

1. Determine stability for: (a) $x' = -x,\ y' = -2y$ ; (b) $x' = y,\ y' = -\sin x - 0.5y$ .
2. Use $V = x^2 + y^2$ to analyze $x' = -x + y^2,\ y' = -y$ .
3. Find all bifurcation points of $\dot x = rx - x^3$ .

Advanced.

4. Prove total energy is a Lyapunov function for the damped pendulum, then apply LaSalle.
5. Analyze the Van der Pol oscillator: show the origin is unstable but a stable limit cycle exists.
6. For Lotka-Volterra competition, derive the conditions for coexistence vs. exclusion.

Programming.

7. Auto-classify 2D linear-system equilibria from trace and determinant; reproduce fig 3.
8. Animate the Hopf bifurcation: sweep $\mu$ and watch the limit cycle grow.

Numerical Side: Letting SciPy Solve the Lyapunov Equation#

On paper, $A^\top P + PA = -Q$ looks like a matrix equation you can stare at. By hand, two dimensions is the limit. Past that, call scipy.linalg.solve_continuous_lyapunov:

1
2
3
4
5
6
7
8

import numpy as np
from scipy.linalg import solve_continuous_lyapunov, eigvals

A = np.array([[0, 1], [-2, -3]])
Q = np.eye(2)
P = solve_continuous_lyapunov(A, -Q)   # convention: A P + P A^T = -Q
print("P =", P)
print("Positive definite?", np.all(eigvals(P).real > 0))

Under the hood it runs Bartels-Stewart: Schur-decompose $A = U T U^\top$ , transform the unknown to $\tilde P = U^\top P U$ , and back-substitute column by column on the resulting upper-triangular system. Cost is $O(n^3)$ , which beats the naive Kronecker-product flattening by a factor of $n^3$ .

Failure mode. If $A$ has eigenvalues on the imaginary axis the equation is singular and SciPy raises LinAlgError. That error is the numerical fingerprint of a non-hyperbolic equilibrium where Hartman-Grobman fails — same phenomenon, two languages.

When Linearization Fails: Center Manifolds#

Hartman-Grobman needs every Jacobian eigenvalue to have non-zero real part. The moment one sits on the imaginary axis, the linear part only tells you the neutral direction; nonlinear terms decide stability.

Toy case: $\dot x = -x^3$ . The Jacobian at $0$ is $0$ , so linear analysis is silent. But $V(x) = x^2/2$ gives $\dot V = -x^4 \le 0$ , hence stable. Lyapunov sees what eigenvalues miss.

Centre Manifold Theorem. With $n_s$ stable and $n_c$ center directions, an invariant $n_c$ -dimensional manifold $W^c$ exists; nearby trajectories collapse onto it exponentially. Stability of the full system reduces to stability on $W^c$ . Hopf normal forms (next chapter) are derived this way.

Practical rule. If np.linalg.eigvals(J) returns an eigenvalue with abs(real) < 1e-8, do not declare stability from spectrum alone. Either build a Lyapunov function or expand to higher Taylor order.

ML Connection: Lyapunov Neural Networks and Stability-Constrained Training#

Two implementation tricks I have actually used:

1. Bake positivity into the architecture. Write $V_\theta(x) = \|\phi_\theta(x)\|^2 + \epsilon \|x\|^2$ with $\phi_\theta$ a standard MLP. Then $V_\theta \ge \epsilon\|x\|^2$ holds by construction, no penalty term needed.
2. Verify the descent condition with an SMT solver. Sample-based loss is not a proof. After training, hand $V_\theta$ and $f$ to dReal or Z3, get a counterexample, retrain — CEGIS loop.

References worth tracking: Neural Lyapunov Control (Chang et al., 2019) and Learning Stability Certificates from Data (Boffi et al., 2020). Already deployed on legged robot controllers.

Summary#

References#

• Strogatz, Nonlinear Dynamics and Chaos, CRC Press (2015)
• Khalil, Nonlinear Systems, Prentice Hall (2002)
• Guckenheimer & Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations, Springer (1983)
• Perko, Differential Equations and Dynamical Systems, Springer (2001)

This is Part 7 of the 18-part series on Ordinary Differential Equations.

• Part 1: Origins and Intuition
• Part 2: First-Order Methods
• Part 3: Higher-Order Linear Theory
• Part 4: The Laplace Transform
• Part 5: Power Series and Special Functions
• Part 6: Linear Systems and the Matrix Exponential
• Part 7: Stability Theory (current)
• Part 8: Nonlinear Systems and Phase Portraits
• Part 9: Chaos Theory and the Lorenz System
• Part 10: Bifurcation Theory
• Part 11: Numerical Methods
• Part 12: Boundary Value Problems
• Part 13: Introduction to PDEs
• Part 14: Epidemic Models
• Part 15: Population Dynamics
• Part 16: Fundamentals of Control Theory
• Part 17: Physics and Engineering Applications
• Part 18: Frontiers and Series Finale