Nonlinear Systems and Linearization

Learning Objectives

By the end of this section you will be able to:

1. Find every fixed point of a 2-D nonlinear system $\mathbf{x}' = \mathbf{F}(\mathbf{x})$ by solving $\mathbf{F}(\mathbf{x}^\ast)=\mathbf{0}$.
2. Compute the Jacobian matrix $J(\mathbf{x}^\ast) = \partial\mathbf{F}/\partial\mathbf{x}\big|_{\mathbf{x}^\ast}$ analytically and verify it numerically.
3. Classify each fixed point as saddle, node, spiral, or center using the trace, determinant, and discriminant of $J$.
4. State and apply the Hartman–Grobman theorem — the rule that says, for hyperbolic fixed points, the nonlinear flow looks topologically like the linearized flow.
5. Identify the failure case (centers, zero eigenvalues) and know that you need a finer tool there.
6. Implement the whole pipeline — fixed-point search, Jacobian, classification — first in NumPy with finite differences, then in PyTorch with exact autograd.

Why Linearize?

“The equations of physics are almost never linear. So why does almost every textbook chapter solve linear systems?”

Because near every interesting point, every smooth nonlinear system looks linear. Zoom in close enough on the orbit of a comet, the swing of a pendulum, the population of a forest, the voltage in a Hodgkin–Huxley neuron, and you will always see the same family of linear pictures we already mastered in sections 23.02–23.06: stable spirals, unstable nodes, saddles, centers.

Linearization is the bridge. It says: solve the linear problem once, then transplant the answer to every fixed point of every smooth nonlinear system in nature. That is the entire content of this section.

The Idea: A Microscope at Each Fixed Point

Picture the flow of a 2-D nonlinear system $\mathbf{x}' = \mathbf{F}(\mathbf{x})$ as a wind field swirling across the plane. At a fixed point $\mathbf{x}^\ast$ the wind dies — air is perfectly still. Drop a leaf at that point, it stays. Drop a leaf slightly off the point, and it drifts. Which way?

Zoom your microscope into a tiny disk around $\mathbf{x}^\ast$. Set $\mathbf{u} = \mathbf{x} - \mathbf{x}^\ast$ (the deviation). A first-order Taylor expansion of $\mathbf{F}$ around $\mathbf{x}^\ast$ says

Drop the higher-order term and the flow $\mathbf{u}' = \mathbf{F}(\mathbf{x}^\ast + \mathbf{u})$ collapses to a linear system:

Step 1 — Finding Fixed Points

A fixed point (a.k.a. equilibrium, critical point, rest point) is any $\mathbf{x}^\ast$ satisfying

In two dimensions this is a system of two equations in two unknowns. Methods:

For the damped pendulum $\dot\theta=\omega,\;\dot\omega=-\sin\theta - b\omega$: $\dot\theta=0$ forces $\omega=0$, then $\dot\omega=0$ with $\omega=0$ gives $\sin\theta=0$, so $\theta\in\{0,\pm\pi,\pm 2\pi,\ldots\}$. The pendulum has infinitely many fixed points — but physically only two are distinct: $(0,0)$ (hanging straight down) and $(\pi,0)$ (balanced straight up).

Step 2 — The Jacobian Matrix

Write the system in components $\mathbf{F}(x,y) = \bigl(F_1(x,y),\;F_2(x,y)\bigr)$. The Jacobian at a point $(x,y)$ is the 2×2 matrix

It is the multivariable generalization of the ordinary derivative $f'(x)$. Its columns tell you what happens to the unit vectors $\mathbf{e}_1, \mathbf{e}_2$ under $\mathbf{F}$; equivalently, its action on any tiny displacement $\mathbf{u}$ is the best linear approximation to $\mathbf{F}(\mathbf{x}^\ast+\mathbf{u}) - \mathbf{F}(\mathbf{x}^\ast)$.

Pendulum example

With $F_1 = \omega$ and $F_2 = -\sin\theta - b\,\omega$ the partials are

So the Jacobian is

Notice $J$ depends only on $\theta$ (and the constant $b$). Evaluated at the two physical fixed points:

The sign of the off-diagonal entry flips because $\cos(\pi) = -1$. That single sign change causes the entire qualitative difference between “pendulum hanging down” and “pendulum balanced up”.

Step 3 — The Linearization Theorem

Once we have $J(\mathbf{x}^\ast)$, the small deviation $\mathbf{u} = \mathbf{x}-\mathbf{x}^\ast$ satisfies the linear system

In plain English: the linearization tells the truth about every hyperbolic fixed point. The only loophole is when some eigenvalue has zero real part — pure imaginary (a center) or exactly zero (a degenerate direction). In those cases the missing higher-order terms decide the actual behavior.

Classifying Fixed Points by Trace and Determinant

For a 2×2 matrix $J$ with trace $\tau = \mathrm{tr}\,J$ and determinant $\Delta = \det J$, the characteristic polynomial is

so the eigenvalues are

Three numbers — $\tau,\,\Delta,\;\tau^2-4\Delta$ — completely determine the fixed-point type. The full chart:

Interactive: Linear Phase-Plane Classifier

Before we attack the nonlinear pendulum, build muscle memory for what each linear type LOOKS like. Drag the four sliders to change the entries of $A$; the panel on the right shows which region of the trace-determinant plane you have landed in. Try the preset buttons — saddle, stable spiral, center — and notice how the field flips between in/out, rotation, decay.

The Damped Pendulum, End to End

Now put the three steps together on the model nonlinear system of this section. Take the damped pendulum

Fixed points

From $\dot\theta=0$: $\omega=0$. From $\dot\omega=0$: $\sin\theta=0$. Physical solutions (mod $2\pi$):

Jacobian

Linearization at the “down” fixed point

$J(0,0) = \begin{bmatrix}0 & 1\\ -1 & -0.3\end{bmatrix}.$ Trace $\tau = -0.3$, determinant $\Delta = 1$, discriminant $\tau^2-4\Delta = 0.09 - 4 = -3.91 < 0$. Eigenvalues:

Complex pair, negative real part ⇒ stable spiral. The pendulum, released from any small angle, will oscillate while losing amplitude. Period of one oscillation: $T = 2\pi / 0.989 \approx 6.35$ time units (matches the period $2\pi/\omega_0$ we know from the small-angle approximation).

Linearization at the “up” fixed point

$J(\pi,0) = \begin{bmatrix}0 & 1\\ +1 & -0.3\end{bmatrix}.$ Trace $\tau = -0.3$, determinant $\Delta = -1 < 0$. Δ negative is the signature of a saddle — no further work needed. Eigenvalues:

One positive eigenvalue (unstable direction) and one negative eigenvalue (stable direction). The unstable manifold direction is the eigenvector for $\lambda = +0.861$: $\mathbf{v}_1 \approx (0.758,\,0.653)$. The stable manifold direction: $\mathbf{v}_2 \approx (-0.653,\,0.758)$. These are the two dashed lines you will see in the interactive picture below.

Glue the two local pictures together: orbits spiral into the down position, except those landing on the stable manifold of the up position, which slide in along $\mathbf{v}_2$; any orbit slightly off the stable manifold gets ejected along $\mathbf{v}_1$ and eventually spirals into the next down position one full revolution away. That entire global story comes from two 2×2 matrices.

Interactive: Pendulum vs Its Linearization

Below: the full nonlinear pendulum field as faint blue arrows; the Jacobian linearization around the highlighted fixed point as bright amber arrows inside a dashed circle. Click anywhere to drop a starting state. The solid bright curve is the TRUE pendulum trajectory; the dashed curve is what the LINEARIZATION predicts from the same starting state. Near the chosen fixed point the two are indistinguishable. Far away — especially across $\theta = \pi/2$ — they diverge dramatically.

Worked Example (Step-by-Step)

Let's do every step by hand, exactly as the computer would, for the system in the simulator ($b = 0.3$).

Computation in Python

Two helpers do all the work: a numerical Jacobian by central differences, and a classifier that reads three numbers off any 2×2 matrix. We loop over the three fixed points of the pendulum and print a one-line verdict for each.

1import numpy as np
2
3# ----------------------------------------------------------------------
4# Damped pendulum:  theta' = omega,  omega' = -sin(theta) - b * omega
5# ----------------------------------------------------------------------
6def F(state, b):
7    theta, omega = state
8    return np.array([omega, -np.sin(theta) - b * omega])
9
10def find_fixed_points():
11    """Return all (theta*, omega*) where F = 0 inside [-pi, pi]."""
12    # omega* = 0  forces  sin(theta*) = 0   ->   theta* in {0, pi, -pi}
13    return [(0.0, 0.0), (np.pi, 0.0), (-np.pi, 0.0)]
14
15def jacobian(state, b, eps=1e-6):
16    """Numerical Jacobian by central differences."""
17    n = len(state)
18    J = np.zeros((n, n))
19    for j in range(n):
20        e = np.zeros(n); e[j] = eps
21        J[:, j] = (F(state + e, b) - F(state - e, b)) / (2 * eps)
22    return J
23
24def classify(J):
25    tr  = np.trace(J)
26    det = np.linalg.det(J)
27    disc = tr**2 - 4*det
28    eigs = np.linalg.eigvals(J)
29    if disc > 0:
30        if det < 0:        name = "saddle"
31        elif tr < 0:       name = "stable node"
32        else:              name = "unstable node"
33    elif disc < 0:
34        if abs(tr) < 1e-9: name = "center (linearization fails)"
35        elif tr < 0:       name = "stable spiral"
36        else:              name = "unstable spiral"
37    else:                  name = "degenerate"
38    return name, eigs
39
40b = 0.3
41for fp in find_fixed_points():
42    J = jacobian(np.array(fp), b)
43    name, eigs = classify(J)
44    print(f"fp = ({fp[0]:+.4f}, {fp[1]:+.4f})  ->  {name}")
45    print(f"   J = {J.tolist()}")
46    print(f"   eigenvalues = {np.array2string(eigs, precision=4)}")

Expected output:

fp = (+0.0000, +0.0000)  ->  stable spiral
   J = [[0.0, 1.0], [-1.0, -0.3]]
   eigenvalues = [-0.15+0.9887j -0.15-0.9887j]
fp = (+3.1416, +0.0000)  ->  saddle
   J = [[0.0, 1.0], [1.0, -0.3]]
   eigenvalues = [ 0.8612 -1.1612]
fp = (-3.1416, +0.0000)  ->  saddle
   J = [[0.0, 1.0], [1.0, -0.3]]
   eigenvalues = [ 0.8612 -1.1612]

PyTorch: Autograd Jacobian and Stability

Finite differences are great for sanity-checking, but they introduce truncation error and require N+1 function evaluations per Jacobian column. PyTorch's $\text{autograd}$ gives the exact Jacobian with a single backward pass — and, crucially, it keeps every quantity inside a differentiable graph. That means we can later optimize stability itself: pick the damping $b$ that puts the pendulum's eigenvalues exactly where we want them.

1import torch
2
3# Damped pendulum vector field, batched over states.
4def F(state, b=0.3):
5    theta, omega = state[..., 0], state[..., 1]
6    return torch.stack([omega, -torch.sin(theta) - b * omega], dim=-1)
7
8def jacobian_autograd(state, b=0.3):
9    """Exact Jacobian via torch.autograd.functional.jacobian — no eps, no error."""
10    return torch.autograd.functional.jacobian(lambda s: F(s, b), state)
11
12def classify(J):
13    eigs = torch.linalg.eigvals(J)
14    real = eigs.real
15    imag = eigs.imag
16    if torch.all(real < 0):
17        return ("stable spiral" if torch.any(imag.abs() > 1e-9)
18                else "stable node")
19    if torch.all(real > 0):
20        return ("unstable spiral" if torch.any(imag.abs() > 1e-9)
21                else "unstable node")
22    if (real[0] * real[1] < 0).item():
23        return "saddle"
24    return "center / degenerate"
25
26pi  = torch.tensor(torch.pi)
27fps = [
28    torch.tensor([0.0,       0.0]),  # spiral
29    torch.tensor([pi.item(), 0.0]),  # saddle
30]
31
32for fp in fps:
33    J     = jacobian_autograd(fp, b=0.3)
34    eigs  = torch.linalg.eigvals(J)
35    name  = classify(J)
36    print(f"fp = ({fp[0]:+.4f}, {fp[1]:+.4f})  ->  {name}")
37    print(f"   J =\n{J}")
38    print(f"   eigenvalues = {eigs}")

Expected output:

fp = (+0.0000, +0.0000)  ->  stable spiral
   J =
tensor([[ 0.0000,  1.0000],
        [-1.0000, -0.3000]])
   eigenvalues = tensor([-0.1500+0.9887j, -0.1500-0.9887j])
fp = (+3.1416, +0.0000)  ->  saddle
   J =
tensor([[ 0.0000,  1.0000],
        [ 1.0000, -0.3000]])
   eigenvalues = tensor([ 0.8612+0.j, -1.1612+0.j])

When Linearization Fails — The Center Case

Hartman–Grobman only protects hyperbolic fixed points (every eigenvalue has nonzero real part). The undamped pendulum $b = 0$ at the origin has Jacobian

with eigenvalues $\lambda = \pm i$ — pure imaginary. Linearization predicts a center (perfect circular orbits around the origin). What does the true nonlinear pendulum do? It conserves the energy $E = \tfrac{1}{2}\omega^2 + (1-\cos\theta)$, and the level curves of $E$ are not circles — they bulge as $\theta$ approaches $\pm\pi$ and connect through the two saddle points to form the famous pendulum separatrix. So the linearization is qualitatively right (closed orbits) but quantitatively wrong (shape), and it cannot say anything at all about the global heteroclinic loop.

Switch on a tiny $b > 0$ and the eigenvalues instantly pick up a small negative real part — the center becomes a very-slowly decaying spiral. Hartman–Grobman is restored, but the warning stands: at the boundary, higher-order terms decide. Tools you meet later — Lyapunov functions, normal forms, center-manifold theory — exist precisely to plug this gap.

Where Linearization Shows Up

Every entry is the same three steps: write $\mathbf{F}$, find $\mathbf{x}^\ast$, eigenvalue-classify $J(\mathbf{x}^\ast)$. The same 2×2 characteristic polynomial decides whether a quadcopter hovers, a forest persists, or an ice age begins.

Common Pitfalls

Summary