Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Write the general solution of $\mathbf{x}'(t) = A\,\mathbf{x}(t)$ when the 2×2 matrix $A$ has two real distinct eigenvalues.
Explain why the eigenvectors are the "natural axes" of the system and how they decouple the equations.
Recognise the three qualitative pictures — stable node, unstable node, saddle — just from the signs of $\lambda_1, \lambda_2$ .
Fit the constants $c_1, c_2$ to an initial condition $\mathbf{x}(0) = \mathbf{x}_0$ by solving a single 2×2 linear system.
Implement the eigenvalue method in NumPy and verify the closed-form solution with PyTorch autograd.

The Setup: Two Coupled Equations

A planar linear system is just two scalar ODEs that mix into each other:

x_1'(t) = a\,x_1(t) + b\,x_2(t), \qquad x_2'(t) = c\,x_1(t) + d\,x_2(t).

Bundling the two unknowns into the vector $\mathbf{x} = (x_1, x_2)^{\top}$ compresses the whole system into one tidy line:

\mathbf{x}'(t) = A\,\mathbf{x}(t), \qquad A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}.

The coupling is the whole problem. If $b = c = 0$ the two equations are independent and we already know the answer from the scalar case: $x_k(t) = x_k(0)\,e^{\lambda_k t}$ .

For generic $A$ the off-diagonal entries couple $x_1$ and $x_2$ : you cannot integrate one without knowing the other. The whole eigenvalue method exists to un-couple the system — to find a new pair of coordinates in which the equations are independent again.

Where the trial solution comes from

Inspired by the scalar case, we look for solutions of the form $\mathbf{x}(t) = e^{\lambda t} \mathbf{v}$ for some scalar $\lambda$ and some constant vector $\mathbf{v}$ . Plug into $\mathbf{x}' = A\mathbf{x}$ :

\lambda\,e^{\lambda t}\,\mathbf{v} = A\,e^{\lambda t}\,\mathbf{v}.

Cancel $e^{\lambda t}$ (never zero) and you get the eigenvalue equation: $A\mathbf{v} = \lambda\mathbf{v}$ .

The Key Insight: Eigenvectors Decouple the System

Suppose $A$ has two real distinct eigenvalues $\lambda_1 \neq \lambda_2$ with eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$ . The trick is to look at the problem in the basis $\{\mathbf{v}_1, \mathbf{v}_2\}$ instead of $\{\mathbf{e}_1, \mathbf{e}_2\}$ .

Any vector can be written as a unique combination $\mathbf{x} = z_1 \mathbf{v}_1 + z_2 \mathbf{v}_2$ . Differentiating in time and using $A\mathbf{v}_k = \lambda_k \mathbf{v}_k$ :

A\mathbf{x} = z_1\,A\mathbf{v}_1 + z_2\,A\mathbf{v}_2 = \lambda_1 z_1\,\mathbf{v}_1 + \lambda_2 z_2\,\mathbf{v}_2,

and matching against $\mathbf{x}' = z_1'\,\mathbf{v}_1 + z_2'\,\mathbf{v}_2$ gives:

z_1'(t) = \lambda_1\,z_1(t), \qquad z_2'(t) = \lambda_2\,z_2(t).

What just happened

In the eigenbasis the two equations are completely independent scalar ODEs. Each one is just exponential growth or decay at its own rate $\lambda_k$ :

z_k(t) = c_k\,e^{\lambda_k t}, \qquad c_k = z_k(0).

Translating back into the original coordinates, $\mathbf{x}(t) = z_1(t) \mathbf{v}_1 + z_2(t) \mathbf{v}_2$ .

The General Solution Formula

Putting it together, every solution of the coupled system is:

\mathbf{x}(t) = c_1\,e^{\lambda_1 t}\,\mathbf{v}_1 \;+\; c_2\,e^{\lambda_2 t}\,\mathbf{v}_2.

Two arbitrary constants $c_1, c_2$ — exactly what we expect for a second-order problem (two scalar equations, two free constants). They are pinned by the initial condition $\mathbf{x}(0) = \mathbf{x}_0$ , which gives the linear system:

\begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} = \mathbf{x}_0, \quad\text{equivalently}\quad V\,\mathbf{c} = \mathbf{x}_0.

Because the two eigenvectors are linearly independent (this is exactly where the "distinct eigenvalues" hypothesis bites), the matrix $V = [\,\mathbf{v}_1\;\mathbf{v}_2\,]$ is invertible and $\mathbf{c} = V^{-1}\mathbf{x}_0$ is unique.

The recipe in five steps

Compute the characteristic polynomial $\det(A - \lambda I) = \lambda^2 - (\operatorname{tr} A)\lambda + \det A = 0$ .
Find the two roots $\lambda_1, \lambda_2$ . Confirm they are real and distinct $(\operatorname{tr} A)^2 - 4\det A > 0$ .
For each $\lambda_k$ , solve $(A - \lambda_k I)\mathbf{v}_k = \mathbf{0}$ to get an eigenvector $\mathbf{v}_k$ .
Write $\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2$ .
Solve the 2×2 linear system $V\mathbf{c} = \mathbf{x}_0$ for the constants.

Three Personalities: Stable Node, Unstable Node, Saddle

The qualitative behavior of every real-distinct-eigenvalue system is controlled by just the signs of $\lambda_1$ and $\lambda_2$ . Three cases, three pictures:

Case	Sign pattern	Origin behavior	Long-term direction
Stable node	λ₁ < λ₂ < 0	All trajectories spiral into 0	Tangent to the SLOWER eigenvector (λ₂, less negative)
Unstable node	0 < λ₁ < λ₂	All trajectories fly out from 0	Tangent to the FASTER eigenvector (λ₂, more positive)
Saddle	λ₁ < 0 < λ₂	Some come in, some go out	Ejected along the unstable eigenvector v₂

Why the "slow direction" wins (nodes)

In a stable node the trajectory is $\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2$ . If $\lambda_1 < \lambda_2 < 0$ , then $e^{\lambda_1 t}$ decays faster than $e^{\lambda_2 t}$ . As $t \to +\infty$ , the $\mathbf{v}_2$ term hangs on longer, so $\mathbf{x}(t) \approx c_2 e^{\lambda_2 t} \mathbf{v}_2$ . The trajectory approaches the origin tangent to the slower eigenvector.

Same algebra in reverse for the unstable node, with $t \to +\infty$ dominated by the faster exponential.

Why the saddle is so dramatic

With $\lambda_1 < 0 < \lambda_2$ , one mode decays and the other grows. Any initial condition with even the tiniest $c_2 \neq 0$ component eventually shoots off along $\mathbf{v}_2$ — because the growing exponential wins, no matter how small $c_2$ starts.

The only way to actually approach the origin is to start exactly on the stable eigenvector line $\mathbf{v}_1$ , where $c_2 = 0$ . That line is called the stable manifold; the opposite line $\mathbf{v}_2$ is the unstable manifold. Saddle points are unforgiving: they are the reason "balancing a pencil on its tip" is impossible.

A handy classifier (no eigenvalues needed)

From the trace $T = \operatorname{tr} A$ and the determinant $D = \det A$ you can read off the case directly without solving the quadratic:

$D < 0$ : saddle (one eigenvalue positive, one negative).
$D > 0,\; T^2 > 4D,\; T < 0$ : stable node.
$D > 0,\; T^2 > 4D,\; T > 0$ : unstable node.

(For the other cases — complex or repeated eigenvalues — see sections 23.05 and 23.06.)

Worked Example: A Saddle Point

Solve

\mathbf{x}'(t) = \begin{bmatrix} 1 & 2 \\ 3 & 0 \end{bmatrix}\mathbf{x}(t), \qquad \mathbf{x}(0) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}.

Click to expand the full hand calculation

Step 1 — Characteristic polynomial.

$A - \lambda I = \begin{bmatrix} 1-\lambda & 2 \\ 3 & -\lambda \end{bmatrix}$ , so

\det(A - \lambda I) = (1-\lambda)(-\lambda) - 6 = \lambda^2 - \lambda - 6.

Cross-check by trace–determinant: $T = 1 + 0 = 1$ , $D = 1\cdot 0 - 2\cdot 3 = -6$ , so $\lambda^2 - T\lambda + D = \lambda^2 - \lambda - 6$ ✓.

Step 2 — Solve for the eigenvalues.

Discriminant: $\Delta = 1 + 24 = 25$ , so $\lambda = (1 \pm 5)/2$ :

\lambda_1 = 3, \qquad \lambda_2 = -2.

Opposite signs ⇒ saddle.

Step 3 — Eigenvector for $\lambda_1 = 3$ .

Solve $(A - 3I)\mathbf{v} = \begin{bmatrix} -2 & 2 \\ 3 & -3 \end{bmatrix}\mathbf{v} = \mathbf{0}$ . Both rows reduce to $-v_1 + v_2 = 0$ , so any nonzero multiple of $(1, 1)$ works.

\mathbf{v}_1 = (1, 1)^{\top}.

Step 4 — Eigenvector for $\lambda_2 = -2$ .

$(A + 2I)\mathbf{v} = \begin{bmatrix} 3 & 2 \\ 3 & 2 \end{bmatrix}\mathbf{v} = \mathbf{0}$ . Both rows give $3v_1 + 2v_2 = 0$ , so $(2, -3)$ (and any scalar multiple) works.

\mathbf{v}_2 = (2, -3)^{\top}.

Step 5 — Write the general solution.

\mathbf{x}(t) = c_1\,e^{3t}\begin{bmatrix} 1 \\ 1 \end{bmatrix} + c_2\,e^{-2t}\begin{bmatrix} 2 \\ -3 \end{bmatrix}.

Step 6 — Apply the initial condition.

At $t = 0$ : $c_1 (1, 1) + c_2 (2, -3) = (1, 0)$ . That gives the 2×2 system $c_1 + 2c_2 = 1$ , $c_1 - 3c_2 = 0$ . Subtract: $5c_2 = 1$ , so $c_2 = \tfrac{1}{5} = 0.2$ and $c_1 = 3c_2 = \tfrac{3}{5} = 0.6$ .

Step 7 — Final closed form.

\boxed{\;\mathbf{x}(t) = 0.6\,e^{3t}\begin{bmatrix} 1 \\ 1 \end{bmatrix} + 0.2\,e^{-2t}\begin{bmatrix} 2 \\ -3 \end{bmatrix}\;}

Step 8 — Sanity-check the IC.

$\mathbf{x}(0) = 0.6(1, 1) + 0.2(2, -3) = (0.6 + 0.4,\; 0.6 - 0.6) = (1, 0)$ ✓.

Step 9 — Watch the saddle take over.

At $t = 0.3$ : $0.6\,e^{0.9} \approx 1.476$ and $0.2\,e^{-0.6} \approx 0.110$ , giving $\mathbf{x}(0.3) \approx (1.695, 1.146)$ .

At $t = 1$ : $0.6\,e^{3} \approx 12.05$ but $0.2\,e^{-2} \approx 0.027$ , so $\mathbf{x}(1) \approx 12.05\,(1, 1) + 0.027\,(2, -3) \approx (12.10, 11.97)$ .

By $t = 1$ the stable mode has shrunk by a factor of $\sim 450$ relative to the unstable one. The trajectory is effectively on the line $\mathbf{v}_1 = (1, 1)$ — exactly the picture the saddle predicts.

Step 10 — Cross-check with the residual.

Differentiate: $\mathbf{x}'(t) = 0.6\cdot 3\,e^{3t}(1, 1) + 0.2\cdot(-2)\,e^{-2t}(2, -3)$ . Compute $A\mathbf{x}(t)$ on each term: $A(1,1) = (3, 3) = 3\,\mathbf{v}_1$ and $A(2,-3) = (-4, 6) = -2\,\mathbf{v}_2$ ✓. Mode by mode the equation holds — and the PyTorch residual below confirms it to machine precision.

Interactive Mode Decomposition

The next playground makes the formula $\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2$ visible. The left panel is the phase plane with the two eigenvector axes (red and blue dashes). The right panel draws the two modes as tip-to-tail arrows and sums them into the purple $\mathbf{x}(t)$ . The bottom strip is the time series of the two components.

Loading mode decomposition explorer…

Four things to discover

Load the Saddle preset and press play. Watch the red mode (unstable) explode while the blue mode (stable) shrinks to nothing. By $t \approx 2$ the purple arrow is essentially on the red eigenvector line — that is the unstable manifold taking over.
On the saddle preset, drag the white initial point exactly onto the blue eigenvector line (you should see $c_1 \approx 0$ in the readout). Press play: the trajectory now actually approaches the origin — that is the stable manifold. Bump the IC off this line by even a hair: the unstable mode wakes up and the trajectory eventually shoots away.
Load Stable node. Both eigenvalues are negative; every trajectory comes in. As $t$ grows the trajectory aligns with the eigenvector for the less-negative eigenvalue — the slow direction wins, the fast direction dies first. Compare the time-series slopes: the red mode (more negative) decays steeper.
Slide the angle between $\mathbf{v}_1$ and $\mathbf{v}_2$ down toward $30^\circ$ . The eigenvectors become nearly parallel, and the matrix $V = [\mathbf{v}_1\;\mathbf{v}_2]$ becomes nearly singular. Notice how $c_1, c_2$ blow up — small angles between eigenvectors are the numerical pain point of this method.

Phase Portrait Playground

For a broader sandbox, here is the general phase-portrait explorer you also met in section 23.03. Drag the four matrix entries directly; click anywhere to drop a new trajectory. Whenever the matrix has real distinct eigenvalues the two dashed eigenvector lines appear, and the info card reads off the classification.

Loading phase-portrait explorer…

Things to play with

Start with the loaded matrix $A = \begin{bmatrix} 1 & 2 \\ 3 & 0 \end{bmatrix}$ — the same saddle as the worked example. Notice the trajectories initially curve from anywhere in the plane, then all asymptote along the same eigenvector line.
Hit the Stable node preset (a diagonal matrix). The eigenvectors line up with the axes; trajectories are just decaying exponentials in each coordinate.
Slide $d$ from 0 up past $d = 6$ . Watch the determinant cross zero — the moment $\det A$ becomes positive (and trace stays positive), the classification card flips from Saddle to Unstable node. That is a bifurcation in the trace-determinant plane.

Computation in Python

We turn the five-step recipe into a self-contained NumPy function. It eigendecomposes $A$ , projects the initial condition into the eigenbasis, propagates each mode by its own scalar exponential, and reassembles $\mathbf{x}(t)$ . No time-stepping, no integration error — purely linear algebra.

Eigenvalue method, vectorized over the time grid

Closed-form 2×2 solver for the real-distinct case

🐍real_distinct_eig_solver.py

Explanation(10)

Code(45)

1Import NumPy

We need numpy.linalg.eig (eigendecomposition) and numpy.linalg.solve (a 2×2 linear solve to convert the initial condition into eigen-coordinates). No SciPy needed — this is pure linear algebra.

3Function signature

Inputs: a 2×2 matrix A, an initial state x0 ∈ ℝ², and a 1-D time grid t. Output: x(t) sampled on the grid, plus the eigen-data so the caller can inspect it. The method runs in closed form — no time-stepping involved.

13Eigendecompose A

np.linalg.eig returns the eigenvalues (a 1-D array) and a matrix V whose COLUMNS are the eigenvectors. Mathematically we factored A = V Λ V⁻¹ where Λ = diag(λ₁, λ₂). This is the whole game: in the basis of V, the system decouples into two scalar ODEs.

14Guard against complex eigenvalues

If the discriminant of the characteristic polynomial is negative, NumPy returns complex eigenvalues. That is the case for spirals and centers (section-05). We refuse it here so callers do not silently get garbage.

18Guard against repeated eigenvalues

If λ₁ = λ₂, the eigenvectors might span only a 1-D subspace — the matrix V is then singular and V⁻¹ does not exist. That degenerate case needs generalized eigenvectors (section-06).

22Project initial condition into the eigenbasis

We need constants c₁, c₂ such that x(0) = c₁v₁ + c₂v₂. Stacking that as a linear system V c = x0 — and solving — gives c. Numerically np.linalg.solve is more stable than computing V⁻¹ explicitly.

25Build each mode on the entire time grid

Once we have c, the trajectory in eigen-coordinates is just c_k(t) = c_k · exp(λ_k t). NumPy broadcasting puts λ_k along axis 0 (rows) and t along axis 1 (columns); the elementwise multiplication and exp produce a (2, T) array — one row per mode, one column per time sample.

28Map back to original coordinates

x(t) = V · c(t). Matrix-multiplying V (which has shape (2, 2)) against C (shape (2, T)) gives an (2, T) array of x-vectors. We transpose to (T, 2) so each row is one snapshot of x.

36Worked example: a saddle

A = [[1, 2], [3, 0]] has trace 1 and determinant -6, so the discriminant is 1 + 24 = 25. Eigenvalues are (1±5)/2 = {3, -2} — opposite signs, hence a SADDLE. With x₀ = (1, 0) the printout shows the trajectory shooting off along the unstable eigenvector v(λ=3) ≈ (1, 1).

41Inspect the eigen-data

Print the eigenvalues, the V matrix, and the eigen-coordinates c. Notice that |c₁| · e^(λ₁·t_max) is much larger than |c₂| · e^(λ₂·t_max) at t=1 — the unstable mode dominates by orders of magnitude. That asymmetry is exactly why the trajectory hugs v₁ for large t.

35 lines without explanation

1import numpy as np
2
3def solve_linear_system_2d(A, x0, t):
4    """
5    Solve  x'(t) = A x(t),  x(0) = x0  for a 2x2 matrix A with real
6    distinct eigenvalues. Returns x(t) on the time grid as shape (len(t), 2).
7
8    Strategy (eigenvalue method):
9      1) Eigendecompose  A = V diag(lambda) V^{-1}.
10      2) Convert initial condition into eigen-coordinates:  c = V^{-1} x0.
11      3) Each eigen-coordinate evolves independently as  c_k(t) = c_k * exp(lambda_k t).
12      4) Map back:  x(t) = V (c(t)).
13    """
14    eigvals, V = np.linalg.eig(A)            # V columns are eigenvectors
15    if not np.all(np.isreal(eigvals)):
16        raise ValueError("This routine handles only real eigenvalues.")
17    lam = np.real(eigvals)
18    V   = np.real(V)
19    if abs(lam[0] - lam[1]) < 1e-9:
20        raise ValueError("Eigenvalues are repeated; need the section-06 method.")
21
22    # c = V^{-1} x0
23    c = np.linalg.solve(V, x0)
24
25    # Modes on the time grid:  C[k, i] = c_k * exp(lam_k * t_i)
26    C = c[:, None] * np.exp(lam[:, None] * t[None, :])
27
28    # Map each column of C back through V to get x(t)
29    X = V @ C                                # shape (2, len(t))
30    return X.T, lam, V, c                    # (T, 2) for plotting
31
32# Worked example from the section: a SADDLE
33A  = np.array([[1.0, 2.0],
34               [3.0, 0.0]])
35x0 = np.array([1.0, 0.0])
36t  = np.linspace(0, 1.0, 6)
37X, lam, V, c = solve_linear_system_2d(A, x0, t)
38
39print("eigenvalues:", lam)
40print("eigenvectors (columns of V):")
41print(V)
42print("coordinates in eigenbasis  c = V^-1 x0 =", c)
43print()
44for ti, xi in zip(t, X):
45    print(f"t={ti:.2f}   x={xi}")

Why this beats a generic ODE solver here

For a linear system with constant coefficients, a black-box integrator (RK4, Dormand-Prince, …) must take many small time steps and accumulates discretization error. The eigenvalue method gives the exact closed form. The savings become enormous when $|\lambda|$ is large (stiff systems) — the "stiffness penalty" only shows up for time-stepping methods.

PyTorch View: Autograd Residual Check

Hand-derived ODE solutions are easy to get wrong. The cheapest possible test: take your candidate $\mathbf{x}(t)$ , autodifferentiate it to get $\mathbf{x}'(t)$ , and check that the residual $R(t) = \mathbf{x}'(t) - A\mathbf{x}(t)$ is machine-zero on every grid point.

Confirm the closed-form saddle solution with automatic differentiation

🐍system_residual_autograd.py

Explanation(12)

Code(46)

1Import PyTorch

Autograd lets us differentiate any tensor-valued Python function symbolically. That makes residual checking for ODE systems a one-liner per component.

3Higher-order factory

We want one residual builder that captures the matrix A and returns a function residual(t, x_fn) checking any candidate trajectory. Closures are the natural fit. Building in float64 because residuals shrink to ~1e-13; float32 would mask them with round-off.

10Activate autograd on t

We rebuild t as a leaf tensor with requires_grad=True so downstream computations are tracked. Without this, torch.autograd.grad would have nothing to differentiate with respect to.

11Evaluate the candidate

x_fn maps t ∈ ℝ^T into x ∈ ℝ^{T×2}. The two columns are the two components x₁(t), x₂(t) of the system. It could be a closed form, a neural network, or any differentiable Python.

13Differentiate each component

torch.autograd.grad(x[:, k].sum(), t) returns ∂(Σ x_k)/∂t. Because the sum is elementwise in t, this equals dx_k/dt evaluated at each grid point. We do it once per component because the grad sums over the output tensor — taking the gradient of the full vector against t would not give us per-component derivatives.

15Stack into the velocity vector

Now xd has shape (T, 2) — one row per time sample, two columns for (dx₁/dt, dx₂/dt). That mirrors the shape of x exactly, so we can subtract Ax row-by-row.

17Compute Ax for every row

x @ A.T performs (x A^T) row-wise, which equals (A x) for column-vectors x — that is, for each row x_i of x we get A x_i in the corresponding row of the output. Using A.T avoids a per-row Python loop.

18Return the residual

If x_fn really is a solution, every entry should be float64-zero (≤1e-13). If something is wrong — wrong eigenvalue, wrong c-coefficients, wrong sign — the residual tells you WHERE on the t-grid the equation breaks and BY HOW MUCH.

26Define the closed form

Plug the hand-derived constants into the formula x(t) = c₁·e^(λ₁t)·v₁ + c₂·e^(λ₂t)·v₂. Notice we use torch.exp (autodifferentiable) rather than math.exp (Python scalar, not tracked). Stacking the two components gives the (T, 2) layout the residual expects.

35Run the residual on a fine grid

200 collocation points on [0, 1] is more than enough — a smooth closed form would be machine-zero everywhere. The maximum residual prints as ~1e-12 in float64, which is the strongest end-to-end sanity check we can run on the hand-derived formula.

38Read off the initial condition

x(0) = c₁v₁ + c₂v₂ = 0.6(1,1) + 0.2(2,-3) = (1.0, 0.0). Confirms the (c₁, c₂) we computed by hand actually pin x(0) to (1, 0).

40Check the trajectory at t = 1

At t = 1 the unstable mode 0.6·e³ ≈ 12.05 dwarfs the stable mode 0.2·e⁻² ≈ 0.027 by ~450×, so x(1) ≈ 12.05·(1, 1) + 0.027·(2, -3) ≈ (12.10, 11.97). The trajectory is essentially on the v₁ line by t = 1 — exactly the saddle behavior the math predicted.

34 lines without explanation

1import torch
2
3def make_system_residual(A_np):
4    """
5    Build a function  R(t, x_fn) = x'(t) - A x(t)
6    that compares the AUTODIFFED derivative of any candidate x(t) against
7    what the system says it should be. If x_fn really is a solution, R == 0
8    on every t.
9    """
10    A = torch.tensor(A_np, dtype=torch.float64)
11
12    def residual(t, x_fn):
13        t = t.detach().clone().requires_grad_(True)
14        x = x_fn(t)                          # shape (T, 2)
15        # x'_k(t) =  d(x_k)/dt for each column k
16        xd0 = torch.autograd.grad(x[:, 0].sum(), t, create_graph=True)[0]
17        xd1 = torch.autograd.grad(x[:, 1].sum(), t, create_graph=True)[0]
18        xd  = torch.stack([xd0, xd1], dim=1)  # (T, 2)
19        # Ax for every row of x
20        Ax = x @ A.T                          # (T, 2)
21        return xd - Ax
22
23    return residual
24
25# Closed-form solution for the worked saddle example.
26# A = [[1, 2], [3, 0]],  x(0) = (1, 0).
27# Eigenvalues: lambda_1 = 3 (v1 = (1, 1)), lambda_2 = -2 (v2 = (2, -3)).
28# Coordinates: c1 = 0.6, c2 = 0.2.
29def x_closed(t):
30    e1 = 0.6 * torch.exp( 3.0 * t)
31    e2 = 0.2 * torch.exp(-2.0 * t)
32    x0 = e1 * 1.0 + e2 * 2.0
33    x1 = e1 * 1.0 + e2 * (-3.0)
34    return torch.stack([x0, x1], dim=1)       # (T, 2)
35
36R = make_system_residual([[1.0, 2.0], [3.0, 0.0]])
37
38t = torch.linspace(0.0, 1.0, 200, dtype=torch.float64)
39r = R(t, x_closed)
40print("max |residual|        =", r.abs().max().item())
41
42# Read off the IC and confirm
43x_at_0 = x_closed(torch.tensor([0.0], dtype=torch.float64))
44print("x(0) =", x_at_0.squeeze().tolist())   # should be [1.0, 0.0]
45x_at_1 = x_closed(torch.tensor([1.0], dtype=torch.float64))
46print("x(1) =", x_at_1.squeeze().tolist())   # roughly (12.10, 11.97)

Bridge to physics-informed neural networks (PINNs)

The residual you just built is the exact loss that PINNs minimise: they parametrise $\mathbf{x}(t; \theta)$ as a neural network and train $\theta$ so that $R$ is small at sampled collocation points. For linear systems we already have the analytic answer — this section's closed form is the gold standard against which every learned solver is measured.

Common Pitfalls

Applying the IC to one mode only

The initial condition pins both constants at once. Solve the full 2×2 system $V\mathbf{c} = \mathbf{x}_0$ ; do not set $c_1 = (\mathbf{x}_0)_1$ and $c_2 = (\mathbf{x}_0)_2$ unless your eigenvectors happen to be the standard basis (i.e. unless $A$ is already diagonal).

Forgetting that eigenvectors are determined only up to scale

Any nonzero scalar multiple of an eigenvector is again an eigenvector. Different textbooks, different software libraries, and different students will land on different normalizations. The general solution is the same — but the numerical values of $c_1, c_2$ change inversely. Always pair your $\mathbf{c}$ with the specific $\mathbf{v}_k$ you used.

Confusing “node” with “saddle”

A node has both eigenvalues of the same sign (stable if negative, unstable if positive). A saddle has eigenvalues of opposite signs. The single quickest test: look at $\det A$ . $\det A = \lambda_1\lambda_2$ is negative exactly for saddles.

Ignoring near-parallel eigenvectors

If the angle between $\mathbf{v}_1$ and $\mathbf{v}_2$ is small, the change-of-basis matrix $V$ is ill-conditioned and $\mathbf{c} = V^{-1}\mathbf{x}_0$ can be huge with massive opposite-sign cancellation in the formula. This is the same numerical issue you saw on the mode-explorer when you shrank the angle slider — the trajectory is still correct mathematically, but small floating-point perturbations move it a lot. In production, prefer $\texttt{np.linalg.solve}$ over forming $V^{-1}$ explicitly.

Trying this method when the eigenvalues are repeated

If $\lambda_1 = \lambda_2$ , the matrix $A$ may have only one independent eigenvector — you do not get two linearly independent solutions of the form $e^{\lambda t}\mathbf{v}$ . The missing solution has the form $(t\,\mathbf{v} + \mathbf{w})e^{\lambda t}$ with $\mathbf{w}$ a generalized eigenvector. That is section 23.06.

Summary

For a 2×2 linear system $\mathbf{x}' = A\mathbf{x}$ whose matrix has two real distinct eigenvalues:

Every solution is a sum of two independent modes: $\mathbf{x}(t) = c_1 e^{\lambda_1 t}\mathbf{v}_1 + c_2 e^{\lambda_2 t}\mathbf{v}_2$ .
The eigenvectors are the natural axes: in their basis the equations decouple into two scalar exponentials.
The constants $c_1, c_2$ are pinned by solving the single 2×2 system $V\mathbf{c} = \mathbf{x}_0$ .
The qualitative picture is determined by the signs of the eigenvalues: both negative ⇒ stable node; both positive ⇒ unstable node; opposite signs ⇒ saddle.
For long-time behavior the eigenvalue with the largest real part wins. In nodes, that "slow" direction is what trajectories asymptote to; in saddles, the positive eigenvalue's direction is the unstable manifold that swallows almost every trajectory.

Cheat sheet

Quantity	Formula	Geometric meaning
Characteristic polynomial	λ² − (tr A) λ + det A	Roots = eigenvalues.
Discriminant	(tr A)² − 4 det A	Positive ⇒ real distinct eigenvalues (this section).
Eigenvector for λₖ	(A − λₖ I) vₖ = 0	Direction in which A acts by scaling, not rotating.
General solution	x(t) = c₁ e^(λ₁t) v₁ + c₂ e^(λ₂t) v₂	Superposition of two pure exponential modes.
Constants	V c = x₀ with V = [v₁ \| v₂]	Project initial condition into the eigenbasis.
Saddle test	det A < 0	Eigenvalues have opposite signs.
Stable-node test	det A > 0, tr A < 0, (tr A)² > 4 det A	Both eigenvalues real and negative.
Unstable-node test	det A > 0, tr A > 0, (tr A)² > 4 det A	Both eigenvalues real and positive.

The slogan

"In the eigenbasis the system decouples — every trajectory is a tip-to-tail sum of two pure exponential modes, one per eigenvector."

Coming next: Section 23.05 tackles complex eigenvalues. When the discriminant goes negative the exponentials become spirals — the closed form is the same superposition, but read through Euler's identity it turns into damped (or growing) oscillation.