Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognize when a 2×2 linear system $\mathbf{x}' = A\mathbf{x}$ has complex eigenvalues directly from the discriminant of its characteristic polynomial.
Convert a complex solution $e^{(\alpha+i\beta)t}(\mathbf{p}+i\mathbf{q})$ into two independent real solutions using Euler's formula.
Read off the geometry of the phase portrait — stable spiral, unstable spiral, or center — from the sign of $\alpha = \mathrm{Re}(\lambda)$ alone.
Compute the period $T = 2\pi/|\beta|$ of the rotation directly from the imaginary part.
Solve the canonical example $\mathbf{x}' = \begin{bmatrix}1 & -2\\ 2 & 1\end{bmatrix}\mathbf{x}$ from scratch and verify it three ways: by hand, by RK4 in plain Python, and by $e^{At}$ in PyTorch.
Interpret the same machinery as the mathematics of mechanical vibration, RLC circuits, and any other system that oscillates while decaying.

The Big Picture: Oscillation Hidden in Algebra

In section 23.04 we saw what happens when a 2×2 matrix $A$ has two real, distinct eigenvalues. The eigenvectors point along two real lines, and every trajectory is a mixture of two pure exponentials $e^{\lambda_1 t}$ and $e^{\lambda_2 t}$ . Phase portraits look like saddles, nodes, or stars.

But many physical systems don't flow along straight lines — they rotate. A pendulum swings back and forth. A capacitor and inductor swap energy. A planet circles its star. Algebraically, the only way for a system of the form $\mathbf{x}' = A\mathbf{x}$ to rotate is for its eigenvalues to have a nonzero imaginary part. That is the signal that the algebra is hiding a sine and a cosine.

The motivating problem

Consider $\mathbf{x}' = \begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix}\mathbf{x}$ . The characteristic equation $\lambda^2 + 1 = 0$ has solutions $\lambda = \pm i$ — pure imaginary, no real part at all.

If we plug $e^{it}\mathbf{v}$ in mechanically, we get a complex trajectory. But the matrix $A$ is real, and any physical $\mathbf{x}(t)$ we measure has to be real. How do we get a real answer out of complex algebra?

The answer is Euler's formula. Once we extract the real and imaginary parts of the complex solution, we get exactly two independent real solutions: one that looks like $\cos(\beta t)$ times one vector minus $\sin(\beta t)$ times another, and one that does the same with the roles of sine and cosine swapped. The geometric picture is a rotation in the plane, optionally scaled by $e^{\alpha t}$ to make it spiral.

When Do Complex Eigenvalues Appear?

For a 2×2 real matrix $A = \begin{bmatrix}a & b\\ c & d\end{bmatrix}$ the characteristic polynomial is

\lambda^2 - \tau\,\lambda + \Delta = 0, \qquad \tau = a+d,\quad \Delta = ad - bc.

The quadratic formula gives

\lambda = \dfrac{\tau \pm \sqrt{\tau^2 - 4\Delta}}{2}.

The eigenvalues are complex precisely when the discriminant $\tau^2 - 4\Delta$ is negative. In that case

\lambda = \underbrace{\dfrac{\tau}{2}}_{=\alpha} \;\pm\; i\,\underbrace{\dfrac{\sqrt{4\Delta - \tau^2}}{2}}_{=\beta}.

The trace–determinant test in one line

Complex eigenvalues ⇔ $\tau^2 < 4\Delta$ . Furthermore $\alpha = \tau/2$ and $\beta^2 = \Delta - \alpha^2 > 0$ , so the determinant must be positive whenever the eigenvalues are complex. This is the strict subset of the trace–determinant plane that sits above the parabola $\Delta = \tau^2/4$ .

Two real numbers $\alpha$ and $\beta$ appear naturally: $\alpha = \mathrm{Re}(\lambda)$ controls growth or decay, and $\beta = \mathrm{Im}(\lambda)$ controls how fast the trajectory rotates. They are the two knobs the rest of the section will turn.

The Complex Solution and Why It's a Problem

The eigenvector method works mechanically even when the eigenvalues are complex. If $\lambda = \alpha + i\beta$ and $A\mathbf{v} = \lambda\mathbf{v}$ , then a formal solution of $\mathbf{x}' = A\mathbf{x}$ is

\mathbf{z}(t) = e^{\lambda t}\,\mathbf{v} = e^{(\alpha+i\beta)t}\,\mathbf{v}.

Two pieces of trouble lurk here.

$\mathbf{v}$ is a complex vector. We cannot point at it on a real plane.
$e^{i\beta t}$ is a complex number for every $t \neq 0$ . So $\mathbf{z}(t)$ has imaginary components, but a physical state vector $\mathbf{x}(t)$ (like position and velocity) must be real.

We are saved by two facts that fit together perfectly:

Real-matrix shortcut. Because $A$ is real, eigenvalues come in conjugate pairs — if $\lambda = \alpha + i\beta$ is one, then $\bar\lambda = \alpha - i\beta$ is the other. Their eigenvectors are also conjugates: $\bar{\mathbf{v}}$ .
Real-linear closure. If $\mathbf{z}(t)$ is a complex solution of a real linear ODE, then $\mathrm{Re}\,\mathbf{z}(t)$ and $\mathrm{Im}\,\mathbf{z}(t)$ are each a real solution of that same ODE.

Why real-linear closure holds

A linear operator $L = d/dt - A$ with real coefficients satisfies $L(\mathrm{Re}\,\mathbf{z}) = \mathrm{Re}\,L\mathbf{z}$ and likewise for the imaginary part. So if $L\mathbf{z} = 0$ , both $L\,\mathrm{Re}\,\mathbf{z} = 0$ and $L\,\mathrm{Im}\,\mathbf{z} = 0$ . Splitting a complex solution into its real and imaginary parts gives two real solutions for free.

The Euler Bridge: From Complex to Real

Write the complex eigenvector as $\mathbf{v} = \mathbf{p} + i\mathbf{q}$ with $\mathbf{p}, \mathbf{q} \in \mathbb{R}^2$ — its real and imaginary parts. Then by Euler's formula $e^{i\beta t} = \cos(\beta t) + i\sin(\beta t)$ ,

\mathbf{z}(t) = e^{\alpha t}\,(\cos\beta t + i\sin\beta t)\,(\mathbf{p} + i\mathbf{q}).

Multiply out the second factor carefully:

(\cos\beta t + i\sin\beta t)(\mathbf{p}+i\mathbf{q}) = \big[\cos\beta t\,\mathbf{p} - \sin\beta t\,\mathbf{q}\big] + i\big[\sin\beta t\,\mathbf{p} + \cos\beta t\,\mathbf{q}\big].

So the real and imaginary parts of $\mathbf{z}$ are

\mathbf{x}_1(t) = e^{\alpha t}\big[\cos(\beta t)\,\mathbf{p} - \sin(\beta t)\,\mathbf{q}\big],

\mathbf{x}_2(t) = e^{\alpha t}\big[\sin(\beta t)\,\mathbf{p} + \cos(\beta t)\,\mathbf{q}\big].

These are two real vector-valued functions, and they are linearly independent (a sine wave and a cosine wave with the same frequency can never be proportional). They form a fundamental set of real solutions for our system.

Two Real Solutions from One Complex Pair

The general real solution is therefore

\boxed{\;\mathbf{x}(t) = e^{\alpha t}\Big[c_1\big(\cos\beta t\,\mathbf{p} - \sin\beta t\,\mathbf{q}\big) + c_2\big(\sin\beta t\,\mathbf{p} + \cos\beta t\,\mathbf{q}\big)\Big].\;}

with two real constants $c_1, c_2$ pinned by the initial condition.

A useful reorganization

Group by basis vector instead of by trig function:

\mathbf{x}(t) = e^{\alpha t}\Big[(c_1\cos\beta t + c_2\sin\beta t)\,\mathbf{p} + (-c_1\sin\beta t + c_2\cos\beta t)\,\mathbf{q}\Big].

The coefficient of $\mathbf{p}$ oscillates 90° out of phase with the coefficient of $\mathbf{q}$ . That phase lag is the algebraic source of the rotation we'll see geometrically next.

Pinning the constants is straightforward. At $t = 0$ :

\mathbf{x}(0) = c_1\,\mathbf{p} + c_2\,\mathbf{q}.

So $(c_1, c_2)$ are the coordinates of the initial state $\mathbf{x}(0)$ in the basis $\{\mathbf{p}, \mathbf{q}\}$ . One 2×2 linear solve and we're done.

Geometric Meaning: Rotation + Scaling

Put the formula in its cleanest form. Choose the basis $\{\mathbf{p}, \mathbf{q}\}$ as your axes in the plane. Then the solution looks like

\begin{bmatrix} u(t) \\ w(t) \end{bmatrix} = e^{\alpha t}\,\begin{bmatrix} \cos\beta t & \sin\beta t \\ -\sin\beta t & \cos\beta t \end{bmatrix}\begin{bmatrix} c_1 \\ c_2 \end{bmatrix}.

Inside the bracket sits a standard rotation matrix by angle $\beta t$ . In front sits a scalar $e^{\alpha t}$ that uniformly stretches or shrinks the plane. That is the entire geometric content of complex eigenvalues:

Complex eigenvalues = rotation + uniform scaling. $\beta$ is the angular frequency of the rotation; $\alpha$ is the exponential rate at which the radius grows ( $\alpha > 0$ ) or shrinks ( $\alpha < 0$ ) or stays constant ( $\alpha = 0$ ).

For the canonical matrix $A = \begin{bmatrix}\alpha & -\beta\\ \beta & \alpha\end{bmatrix}$ — the one the interactive explorer below uses — this geometric picture is exact in the standard $(x, y)$ axes, with $\mathbf{p} = (1, 0)$ and $\mathbf{q} = (0, -1)$ . For a general $A$ with complex eigenvalues, the picture is the same up to a linear change of coordinates into the basis $\{\mathbf{p}, \mathbf{q}\}$ .

The Three Regimes: Spirals and Centers

The sign of $\alpha$ selects one of three qualitatively distinct phase portraits.

Sign of α	Name	Geometry	As t → ∞	Energy interpretation
α < 0	Stable spiral (sink)	Trajectories spiral INTO the origin	‖x‖ → 0	Energy is being dissipated (damping)
α = 0	Center	Trajectories trace CLOSED orbits (ellipses)	‖x‖ stays bounded — pure rotation	Energy is conserved exactly
α > 0	Unstable spiral (source)	Trajectories spiral AWAY from the origin	‖x‖ → ∞	Energy is being injected (driving)

Notice that direction of rotation (clockwise vs counterclockwise) is independent of all this — it depends only on the sign of the off-diagonal entry $c$ in the general matrix (positive $c$ → counterclockwise). The growth/decay regime and the rotation direction are independent knobs.

The trace–determinant map

On the $(\tau, \Delta)$ plane, complex eigenvalues live above the parabola $\Delta = \tau^2/4$ . Within that region:

$\tau > 0$ → unstable spirals,
$\tau = 0$ → centers (the positive $\Delta$ axis),
$\tau < 0$ → stable spirals.

Section 23.03 builds the full map. Complex eigenvalues are the "upper half" of it.

Worked Example: $\mathbf{x}' = \begin{bmatrix}1 & -2\\ 2 & 1\end{bmatrix}\mathbf{x}$

Solve with initial condition $\mathbf{x}(0) = (1, 0)$ .

Click to expand the full hand calculation

Step 1 — Characteristic polynomial.

$A - \lambda I = \begin{bmatrix}1-\lambda & -2\\ 2 & 1-\lambda\end{bmatrix}$ . Determinant $(1-\lambda)^2 + 4 = 0$ , so $(1-\lambda)^2 = -4$ and $\lambda = 1 \pm 2i$ .

So $\alpha = 1$ (growth rate) and $\beta = 2$ (angular frequency). Already we know: unstable spiral, period $T = 2\pi/2 = \pi$ .

Step 2 — Find one complex eigenvector.

Pick $\lambda = 1 + 2i$ . Solve $(A - \lambda I)\mathbf{v} = \mathbf{0}$ :

\begin{bmatrix} -2i & -2 \\ 2 & -2i \end{bmatrix}\begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \mathbf{0}.

The first row gives $-2i\,v_1 - 2\,v_2 = 0$ , so $v_2 = -i\,v_1$ . Choose $v_1 = 1$ to obtain

\mathbf{v} = \begin{bmatrix} 1 \\ -i \end{bmatrix} = \underbrace{\begin{bmatrix}1\\0\end{bmatrix}}_{\mathbf{p}} + i\,\underbrace{\begin{bmatrix}0\\-1\end{bmatrix}}_{\mathbf{q}}.

Quick check on the second row: $2\cdot 1 + (-2i)\cdot(-i) = 2 - 2 = 0$ ✓.

Step 3 — Form the two real solutions.

With $\alpha = 1$ , $\beta = 2$ , $\mathbf{p} = (1, 0)$ , $\mathbf{q} = (0, -1)$ :

\mathbf{x}_1(t) = e^{t}\big[\cos 2t\,(1,0) - \sin 2t\,(0,-1)\big] = e^{t}\,(\cos 2t,\;\sin 2t),

\mathbf{x}_2(t) = e^{t}\big[\sin 2t\,(1,0) + \cos 2t\,(0,-1)\big] = e^{t}\,(\sin 2t,\;-\cos 2t).

Step 4 — General solution.

\mathbf{x}(t) = e^{t}\Big( c_1(\cos 2t,\, \sin 2t) + c_2(\sin 2t,\, -\cos 2t)\Big).

Step 5 — Apply $\mathbf{x}(0) = (1, 0)$ .

At $t = 0$ : $\mathbf{x}(0) = c_1(1, 0) + c_2(0, -1) = (c_1,\, -c_2)$ . Match to $(1, 0)$ : $c_1 = 1$ , $c_2 = 0$ .

Step 6 — Final answer.

\boxed{\;\mathbf{x}(t) = e^{t}\,(\cos 2t,\; \sin 2t).\;}

This describes a point starting at $(1, 0)$ , rotating counterclockwise with angular frequency 2, while its distance from the origin grows like $e^t$ . After one period $t = \pi$ it has returned to angle 0 but at radius $e^\pi \approx 23.14$ .

Step 7 — Verify by substitution.

$x'(t) = e^t \cos 2t - 2 e^t \sin 2t = x - 2 y$ ✓ (first row of $A\mathbf{x}$ ).

$y'(t) = e^t \sin 2t + 2 e^t \cos 2t = 2 x + y$ ✓ (second row).

Step 8 — Spot check at $t = 0.5$ .

$e^{0.5} \approx 1.6487$ , $\cos 1 \approx 0.5403$ , $\sin 1 \approx 0.8415$ .

So $x(0.5) \approx 1.6487 \cdot 0.5403 \approx 0.8908$ and $y(0.5) \approx 1.6487 \cdot 0.8415 \approx 1.3874$ . The Python and PyTorch computations below should reproduce these numbers to many decimal places.

Interactive Phase-Portrait Explorer

The explorer below fixes the matrix in its canonical form $A = \begin{bmatrix}\alpha & -\beta\\ \beta & \alpha\end{bmatrix}$ so the two sliders are exactly the real and imaginary parts of the eigenvalues. $\mathbf{p} = (1, 0)$ (gold arrow) and $\mathbf{q} = (0, -1)$ (purple arrow) are the real/imaginary parts of the eigenvector. Click anywhere in the plane to drop a new initial condition; its trajectory is integrated forward (color) and backward (gray) by RK4.

Loading complex-eigenvalues phase portrait…

Experiments to run

Push $\alpha$ all the way to $0$ . Every orbit becomes a perfect circle of radius equal to its initial distance from the origin — the center regime. Now nudge $\alpha$ slightly negative — the circles tilt into inward spirals (stable). Slightly positive — outward spirals (unstable). The qualitative picture flips at the SIGN of $\alpha$ , not its size.
With $\alpha = 0$ , change $\beta$ and watch the animated dot speed up. Halving $\beta$ exactly doubles the period — match it to the displayed $T = 2\pi/|\beta|$ .
Set $\alpha = -0.5$ , $\beta = 3$ . The trajectories make many rotations before noticeably shrinking — typical of a lightly damped oscillator. Pendulums and tuning forks live here.
Set $\alpha = -1$ , $\beta = 0.3$ . The decay dominates the rotation — trajectories barely curl before reaching the origin. This is heavy damping.

The Time-Series View

The phase portrait shows the trajectory as a curve in the plane. The time-series view shows each component $x(t)$ and $y(t)$ separately, as oscillating functions of time. Both views describe the same motion — but the time-series view is what an oscilloscope would record, and it makes the role of $\alpha$ (envelope) versus $\beta$ (frequency) crystal clear.

Loading time-series view…

Two views, one motion

The phase portrait suppresses time; the time-series view suppresses geometry. The two views together are the complete picture. The envelope $\pm e^{\alpha t}$ in the time view is the same thing as the radial growth in the phase portrait. The phase lag between $x(t)$ and $y(t)$ is the same thing as the rotation. Trade fluency between them.

Computing the Solution in Python

Two complementary computations. First the closed-form path: extract the complex eigenvalue, split into $\mathbf{p}, \mathbf{q}$ , assemble the formula, pin the constants. Then the numerical path: integrate the ODE directly with RK4 and verify both branches agree to many digits. The two paths are independent — if they meet in the middle, both are right.

Path 1: closed-form via the eigendecomposition

Build the real general solution from one complex eigenpair

🐍complex_eig_closed_form.py

Explanation(12)

Code(43)

1Import NumPy

We need two things from NumPy. (1) np.linalg.eig for eigenvalues and eigenvectors of a general matrix — it returns COMPLEX dtypes when the eigenvalues are complex, which is exactly what we want. (2) Vectorised scalar functions np.exp, np.cos, np.sin so the entire trajectory on a time grid is computed in one expression.

4Build the coefficient matrix A

This is the matrix from the worked example. Its trace is 1 + 1 = 2 and its determinant is 1·1 − (−2)·2 = 5. The discriminant of the characteristic polynomial is tr² − 4 det = 4 − 20 = −16 < 0, which guarantees complex eigenvalues. We do not have to predict the eigenvalues by hand — eig will compute them.

8Compute eigenvalues and eigenvectors with np.linalg.eig

np.linalg.eig solves det(A − λI) = 0 numerically and returns (1) a 1-D array w of length 2 holding the eigenvalues — here w = [1+2j, 1−2j], appearing as a conjugate pair because A is real; and (2) a 2×2 matrix V whose columns are the corresponding eigenvectors v_1, v_2. When eigenvalues are complex, V is complex too. The columns are paired with the eigenvalues by INDEX, not by name.

13Pick one eigenvalue and its eigenvector

Because A is real, the eigenvalues and eigenvectors come in conjugate pairs (λ, v) and (λ̄, v̄). We only need ONE of them — the other adds no new information. We pick the first column. lam ends up as the scalar 1+2j and v is the complex column vector [0 − 0.707j, −0.707 + 0j].

17Split v into its real and imaginary parts: v = p + i q

This is the crucial step that takes us from complex to real. Any complex eigenvector v can be written as v = p + i q where p = Re(v) and q = Im(v) are both REAL 2-vectors. NumPy's .real and .imag attributes pull them out elementwise. For our v: p = [0, −0.707], q = [−0.707, 0]. Notice neither p nor q is itself an eigenvector — together they span the same real 2D plane as v and v̄.

19Print α, β, p, q to make the structure visible

We read α = Re(λ) = 1 (the growth rate) and β = Im(λ) = 2 (the angular frequency). p and q form the basis for the real plane on which the rotation+scaling happens. Numpy's eigenvectors are unit-normalised in the COMPLEX inner product, so p and q may not be unit length individually — that is fine, it just shows up later as a constant overall scale that the initial condition will absorb.

24general_solution — the closed-form formula

Every real solution of x' = A x with complex eigenvalues α ± iβ has the form e^{αt} times a linear combination of two basis trajectories: cos(βt)·p − sin(βt)·q and sin(βt)·p + cos(βt)·q. These come straight from the real and imaginary parts of the complex solution e^{(α+iβ)t}(p+iq), expanded with Euler's formula. The output is a (T, 2) array — one row per time, two columns for x and y.

28Why the time-grid broadcasting trick (env[:, None])

p and q are shape (2,). cos(βt) and sin(βt) on a time grid are shape (T,). We need to multiply each time sample by the same 2-vector. The trick env[:, None] reshapes (T,) → (T, 1) so that NumPy broadcasts against (2,) into (T, 2). Without the [:, None] the multiplication would fail with a shape-mismatch error or, worse, produce a scalar dot product silently. This is the standard NumPy idiom for 'one row per time'.

35Pin the two constants c1, c2 to the initial condition

At t = 0 we have e^{α·0} = 1, cos(0) = 1, sin(0) = 0. So x(0) = c1 · p + c2 · q. To match the initial condition x(0) = (1, 0), we solve the 2×2 linear system [p | q] · [c1; c2] = (1, 0). M = np.column_stack([p, q]) builds the matrix whose columns are p and q. np.linalg.solve does the rest.

39Use the worked-example initial condition

x(0) = (1, 0) is the IC used in the hand calculation. After this line c1 ≈ 0 and c2 ≈ −1.414 — the exact values depend on how NumPy normalised p and q, but the product c1·p + c2·q reproduces (1, 0) by construction.

41Evaluate the closed-form on a small time grid

t = np.linspace(0, 0.5, 6) is six points from 0 to 0.5. We pass them through general_solution to get the trajectory's value at each sample. The last row xt[-1] is x(0.5).

44Cross-check against the hand-derived answer

From the hand derivation, x(t) = e^t cos(2t), y(t) = e^t sin(2t). At t = 0.5: e^0.5 ≈ 1.6487, cos(1) ≈ 0.5403, sin(1) ≈ 0.8415, so x(0.5) ≈ 0.8908 and y(0.5) ≈ 1.3874. The print should show approximately [0.890808, 1.387351]. If it does, both branches of the derivation agree.

31 lines without explanation

1import numpy as np
2
3# A from the worked example
4A = np.array([[1.0, -2.0],
5              [2.0,  1.0]])
6
7# Step 1: get eigenvalues + eigenvectors (general approach)
8w, V = np.linalg.eig(A)
9print("eigenvalues:", w)        # [1.+2.j  1.-2.j]
10print("eigenvectors:")
11print(V)
12
13# Pick one eigenvalue/eigenvector of the conjugate pair
14lam = w[0]                       # 1 + 2j
15v   = V[:, 0]                    # complex column
16
17# Split v into real + imaginary parts: v = p + i q
18p = v.real
19q = v.imag
20print("alpha =", lam.real, "  beta =", lam.imag)
21print("p =", p, "  q =", q)
22
23# Step 2: build the real general solution
24#   x(t) = e^{alpha t} [ c1 (cos(beta t) p - sin(beta t) q)
25#                      + c2 (sin(beta t) p + cos(beta t) q) ]
26def general_solution(c1, c2, t):
27    a, b = lam.real, lam.imag
28    env  = np.exp(a * t)
29    c, s = np.cos(b * t), np.sin(b * t)
30    # Each row of the result is x(t_i) in R^2
31    return env[:, None] * (c1 * (c[:, None]*p - s[:, None]*q)
32                         + c2 * (s[:, None]*p + c[:, None]*q))
33
34# Step 3: pin (c1, c2) to the initial condition x(0) = (1, 0)
35# At t=0:  x(0) = c1 * p + c2 * q
36M = np.column_stack([p, q])
37c1, c2 = np.linalg.solve(M, np.array([1.0, 0.0]))
38print("c1 =", c1, "  c2 =", c2)
39
40# Step 4: evaluate the answer and cross-check
41t  = np.linspace(0.0, 0.5, 6)
42xt = general_solution(c1, c2, t)
43print("x(0.5) =", xt[-1])

Path 2: pure numerical integration with RK4

Solve the same system without ever computing an eigenvalue

🐍complex_eig_rk4.py

Explanation(11)

Code(35)

3rk4_step — one step of classical Runge–Kutta 4

RK4 evaluates the right-hand side f(x) = A x at FOUR carefully chosen points within the interval [t, t+dt] and averages them with weights (1, 2, 2, 1)/6. The result is a numerical scheme whose error per step is O(dt^5) and whose global error is O(dt^4). For a linear system this is overkill — the exact solution can be written down with eig — but RK4 is the universal tool you reach for once you start studying nonlinear systems later in the chapter.

5k1 — slope at the start of the interval

k1 = A @ x is the time derivative at the current state x. Multiplying by dt would give an Euler step. RK4 uses k1 only as a probe to look further into the interval.

6k2 — slope at the midpoint, using k1 to predict the midpoint state

We tentatively step half-way using k1 (the prediction x + 0.5·dt·k1) and re-evaluate the slope there. This is the corrector idea: don't trust the slope at the start; sample again somewhere closer to where you're going.

7k3 — refined midpoint slope, using k2 to predict

Same idea, but we predict the midpoint using k2 instead of k1. This is the second refinement of the midpoint slope. The two midpoint estimates k2 and k3 are weighted equally in the final average.

8k4 — slope at the end of the interval

Predict x at t + dt using the refined midpoint slope k3, then sample the right-hand side there. k4 is the 'arrival' slope.

9Weighted average — the magic of RK4

x_{new} = x + (dt/6) (k1 + 2 k2 + 2 k3 + k4). The middle two slopes get double weight because they probe the midpoint, which Taylor-series analysis shows is where the most error is hiding. The 1:2:2:1 weighting is exactly what cancels the dt^2, dt^3, and dt^4 error terms.

12simulate — drive RK4 across a time grid

out[i] depends only on out[i-1] and the step size dt — this is the textbook serial structure of any explicit ODE integrator. For a real production solver you would use scipy.integrate.solve_ivp, which adapts dt automatically; here the explicit loop makes the mechanics transparent.

22Reproduce the worked example

Same A and same initial condition (1, 0) as in the hand derivation. 501 grid points on [0, 0.5] gives dt = 0.001 — comfortably small for RK4 on this matrix.

28Build the closed-form trajectory for comparison

From the hand derivation, the exact solution starting at (1, 0) is x(t) = e^t cos(2t), y(t) = e^t sin(2t). np.column_stack glues two (T,) arrays into a (T, 2) array with the same shape as the simulator output. Now we can subtract.

30Max absolute error between numerics and exact

With dt = 0.001, RK4 on this 2×2 problem easily reaches ~10^-12 to 10^-14 max error — limited by floating-point round-off, not by RK4 itself. A first-order Euler scheme on the same grid would give ~10^-3 error, three orders of magnitude worse. This is why we reach for RK4 (or better) for any ODE study.

32Final printout

Both endpoints should agree to many digits with x(0.5) ≈ (0.890808, 1.387351). If you see a mismatch in the first few digits, dt is too large or the matrix entries got typo-edited.

24 lines without explanation

1import numpy as np
2
3def rk4_step(A, x, dt):
4    """One classical RK4 step for the linear system x' = A x."""
5    k1 = A @ x
6    k2 = A @ (x + 0.5 * dt * k1)
7    k3 = A @ (x + 0.5 * dt * k2)
8    k4 = A @ (x + dt * k3)
9    return x + (dt / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
10
11
12def simulate(A, x0, t_grid):
13    """Integrate x' = A x on the given time grid using RK4."""
14    out = np.zeros((len(t_grid), len(x0)))
15    out[0] = x0
16    for i in range(1, len(t_grid)):
17        dt = t_grid[i] - t_grid[i-1]
18        out[i] = rk4_step(A, out[i-1], dt)
19    return out
20
21
22# Reproduce the worked example: x(0) = (1, 0), follow it to t = 0.5
23A  = np.array([[1.0, -2.0],
24               [2.0,  1.0]])
25x0 = np.array([1.0, 0.0])
26t  = np.linspace(0.0, 0.5, 501)
27xt = simulate(A, x0, t)
28
29# Compare against the closed form
30x_exact = np.column_stack([np.exp(t)*np.cos(2*t),
31                            np.exp(t)*np.sin(2*t)])
32err = np.max(np.abs(xt - x_exact))
33print(f"max |x_rk4 - x_exact| = {err:.3e}")
34print("x_rk4(0.5)   =", xt[-1])
35print("x_exact(0.5) =", x_exact[-1])

Why teach both paths

The closed-form path is the structural answer — it tells you what kind of trajectory you have and exposes the two knobs $\alpha, \beta$ . The numerical path is the universal answer — it doesn't care whether the system is linear or whether the eigenvalues even exist. When you move to nonlinear systems in section 23.07 the closed-form path disappears, but the RK4 path still works.

PyTorch View: Eigendecomposition and Matrix Exponential

PyTorch gives us a third independent route to the same answer: $\mathbf{x}(t) = e^{At}\mathbf{x}_0$ , the matrix exponential. This is the most elegant closed form possible — one formula, no special cases, valid whether the eigenvalues are real, complex, or even repeated. We then use autograd to verify that the closed-form trajectory satisfies the ODE pointwise.

Matrix exponential + autograd as the universal ODE verifier

🐍complex_eig_pytorch.py

Explanation(11)

Code(40)

1Import PyTorch

We use three pieces of torch here: (1) torch.linalg.eig, the modern replacement for the deprecated torch.eig, which returns complex eigenvalues directly; (2) torch.linalg.matrix_exp, a vectorised matrix exponential that handles the eA·t in closed form; (3) torch.autograd.grad to differentiate a closed-form trajectory and verify it satisfies the ODE.

4Build A as a real tensor

torch.tensor with floats gives a real (float32) tensor. We do NOT need to cast to complex — torch.linalg.eig is happy to return complex eigenvalues from a real input matrix.

8torch.linalg.eig — modern eigendecomposition

Returns two complex tensors. eigvals has shape (n,) and contains the eigenvalues. eigvecs is the (n, n) matrix whose columns are the eigenvectors. For our A you'll see eigvals = tensor([1+2j, 1-2j]) and eigvecs containing complex columns. PyTorch's API mirrors NumPy's np.linalg.eig — so the structure of the previous Python example transfers verbatim.

15Pick the time t = 0.5

Match the worked example's evaluation point so we can verify against the hand-derived answer (0.8908, 1.3874).

17torch.linalg.matrix_exp — the operator-theoretic solution

The matrix exponential e^{A·t} is the unique solution operator for x' = A x: starting at x_0, the state at time t is exp(A·t) · x_0. PyTorch's matrix_exp is computed by scaling and squaring with a Padé approximation, so it handles non-diagonalizable matrices and complex eigenvalues uniformly — no special case for the complex regime. This makes it the cleanest one-liner for ANY linear ODE.

18Apply the operator to x0

Matrix-vector multiplication. After this line, x_at_half should be tensor([0.8908, 1.3874]) (up to float32 precision, ~1e-6 relative error). Compare against the NumPy and hand-derived numbers — all three branches of the section now agree.

23Define the closed-form solution as a torch callable

x_of_t takes a scalar t and returns the (x, y) tensor at that time. Using torch.exp/torch.cos/torch.sin (not math.exp etc.) is critical — only torch operations participate in autograd. Mixing math.* and torch.* breaks the computation graph silently and gives wrong gradients.

31Make t a leaf tensor with requires_grad

Autograd only differentiates with respect to tensors that have requires_grad=True. We make a fresh leaf tensor at t = 0.7 — a generic interior point, NOT the easy t = 0. From here, every operation involving t_var is recorded into the autograd graph.

33Differentiate component-by-component

torch.autograd.grad(scalar, t_var) returns d(scalar)/dt as a one-element tuple. We do it once per component because x is a 2-vector and grad expects a scalar output (we'd otherwise need to pass grad_outputs). The result xdot is a (2,) tensor whose entries are dx/dt and dy/dt at t = 0.7. create_graph=True on the first call is not strictly needed here, but it would allow a second derivative if we wanted one.

37Compute A · x(0.7) for comparison

If our closed form really solves x' = A x, then for every t we must have x'(t) = A · x(t). We compute the right-hand side A @ x_val with x_val detached (we don't want autograd to backpropagate into x_val again).

40The residual — the gold-standard sanity check

abs(xdot − A x).max() at t = 0.7 should be on the order of 10^-6 to 10^-7 in float32, limited by round-off. If it is large (say 0.1 or more) something is wrong in the closed form, or torch.cos / torch.exp got accidentally replaced by their math.* counterparts. This is the same residual idea you'll see again in the next chapter when we cover physics-informed neural networks.

29 lines without explanation

1import torch
2
3# Same worked example, computed with PyTorch
4A = torch.tensor([[1.0, -2.0],
5                  [2.0,  1.0]])
6
7# 1. Eigendecomposition of A — works on complex eigenvalues too
8eigvals, eigvecs = torch.linalg.eig(A)
9print("eigvals:", eigvals)        # tensor([1+2j, 1-2j])
10print("eigvecs:")
11print(eigvecs)
12
13# 2. Matrix exponential — the cleanest closed-form formula
14#    For x' = A x with x(0) = x0,  x(t) = expm(A t) @ x0.
15t  = torch.tensor(0.5)
16x0 = torch.tensor([1.0, 0.0])
17expAt = torch.linalg.matrix_exp(A * t)
18x_at_half = expAt @ x0
19print("expm(A*0.5) =\n", expAt)
20print("x(0.5)      =", x_at_half)
21
22# 3. Cross-check via autograd: verify x' = A x at one point
23def x_of_t(t_scalar: torch.Tensor) -> torch.Tensor:
24    """Closed-form: x(t) = (e^t cos 2t, e^t sin 2t)."""
25    return torch.stack([
26        torch.exp(t_scalar) * torch.cos(2 * t_scalar),
27        torch.exp(t_scalar) * torch.sin(2 * t_scalar),
28    ])
29
30t_var = torch.tensor(0.7, requires_grad=True)
31x_val = x_of_t(t_var)
32# Per-component derivatives dx_i/dt
33xdot = torch.stack([
34    torch.autograd.grad(x_val[0], t_var, create_graph=True)[0],
35    torch.autograd.grad(x_val[1], t_var)[0],
36])
37ax = A @ x_val.detach()
38print("x'(0.7) via autograd:", xdot.detach())
39print("A x(0.7)             :", ax)
40print("residual             :", (xdot.detach() - ax).abs().max().item())

Why the matrix exponential is the "right" answer

For a linear system $\mathbf{x}' = A\mathbf{x}$ the operator $e^{At}$ propagates ANY initial condition to the state at time $t$ . It packages the eigenvalue analysis (real, complex, repeated) into a single object. When $A$ has complex eigenvalues $\alpha \pm i\beta$ , the matrix $e^{At}$ contains entries built from $e^{\alpha t}\cos\beta t$ and $e^{\alpha t}\sin\beta t$ — exactly the building blocks we derived by hand. Chapter 24 (Laplace transforms) gives you a third way to compute it.

Physical Meaning: Where This Appears in Nature

System	Equations	α and β	Regime
Mass on a spring with damping	m ẍ + c ẋ + k x = 0	α = −c / 2m, β = √(k/m − α²)	Stable spiral (under-damped)
Pendulum (small angle, no friction)	θ̈ + (g/L) θ = 0	α = 0, β = √(g/L)	Center
RLC circuit (series)	L Q̈ + R Q̇ + Q/C = 0	α = −R / 2L, β = √(1/LC − α²)	Stable spiral
Predator-prey near coexistence	linearization of Lotka–Volterra	α = 0 (at the fixed point), β = √(r·s)	Center (locally)
Coupled neurons (simple model)	v̇ = … ẇ = …	depends on synaptic strengths	Stable or unstable spiral

Every one of these is a 2D linear (or linearized) system whose matrix has complex eigenvalues. Whenever you see damped oscillation in the physical world, the underlying algebra is a spiral in some phase plane, and the math you just learned is the complete description of it.

The under-damped oscillator in one line

For $m\ddot x + c\dot x + k x = 0$ with $c^2 < 4mk$ (under-damped), let $\omega_0 = \sqrt{k/m}$ be the undamped frequency and $\zeta = c/(2\sqrt{mk})$ the damping ratio. Then $\alpha = -\zeta\omega_0$ and $\beta = \omega_0\sqrt{1-\zeta^2}$ . The decay envelope is $e^{-\zeta\omega_0 t}$ and the oscillation has period $T = 2\pi / (\omega_0\sqrt{1-\zeta^2})$ . Every formula in this section specializes to this one engineering recipe.

Common Pitfalls

Confusing v with p and q

$\mathbf{v} = \mathbf{p} + i\mathbf{q}$ is the complex eigenvector. Individually, neither $\mathbf{p}$ nor $\mathbf{q}$ is an eigenvector of $A$ . They are real vectors that span the real 2D plane on which the rotation happens. Don't try to verify $A\mathbf{p} = \alpha\mathbf{p}$ — you'll find it doesn't hold.

Forgetting the minus sign in cos(βt) p − sin(βt) q

The real part of $(\cos+i\sin)(\mathbf{p}+i\mathbf{q})$ is $\cos\beta t\,\mathbf{p} - \sin\beta t\,\mathbf{q}$ , with a minus. It comes from $i \cdot i = -1$ . Forgetting the minus gives a solution that satisfies a related but different ODE — the verification step in Step 7 of the worked example will catch it.

Mistaking a center for a spiral on a finite-time plot

On a numerical phase portrait integrated for only a few periods, a weak spiral (small $|\alpha|$ ) and a true center look almost identical. To distinguish, run the simulation for many periods or — better — compute the eigenvalues symbolically. Numerical drift in an RK4 simulation can even make a mathematical center look like a slow spiral; that's a discretization artifact, not physics.

Picking the "wrong" eigenvalue of the pair

If you happen to pick $\bar\lambda = \alpha - i\beta$ with eigenvector $\bar{\mathbf{v}} = \mathbf{p} - i\mathbf{q}$ instead, you get the same two real solutions back — just with $\mathbf{q}$ negated and a sign flip in the formula. The general real solution is the same set. The conjugate eigenpair carries no new information.

Numerical eigenvectors are normalised in the COMPLEX inner product

When you call np.linalg.eig the returned eigenvectors satisfy $\mathbf{v}^* \mathbf{v} = 1$ using the conjugate transpose. So $\|\mathbf{p}\|^2 + \|\mathbf{q}\|^2 = 1$ but neither $\mathbf{p}$ nor $\mathbf{q}$ is individually unit length. That is fine — the overall scaling is absorbed when you fit the constants $c_1, c_2$ to the initial condition.

Summary

Complex eigenvalues are not an exotic edge case — they are the algebra of every oscillating linear system. The recipe is short:

Compute eigenvalues $\lambda = \alpha \pm i\beta$ of $A$ .
Find one complex eigenvector $\mathbf{v} = \mathbf{p} + i\mathbf{q}$ .
Form the real general solution $\mathbf{x}(t) = e^{\alpha t}[c_1(\cos\beta t\,\mathbf{p} - \sin\beta t\,\mathbf{q}) + c_2(\sin\beta t\,\mathbf{p} + \cos\beta t\,\mathbf{q})]$ .
Pin $c_1, c_2$ by solving the 2×2 system $[\mathbf{p}\;\mathbf{q}](c_1, c_2)^\top = \mathbf{x}(0)$ .

Master picture

Geometrically, the flow is a rotation by angle $\beta t$ combined with a uniform radial scaling by $e^{\alpha t}$ . The sign of $\alpha$ alone selects:

$\alpha < 0$ : stable spiral — orbits decay into the origin.
$\alpha = 0$ : center — orbits are closed ellipses.
$\alpha > 0$ : unstable spiral — orbits grow without bound.

The slogan

"Complex eigenvalues are α + iβ; the real solution is a rotation by βt scaled by e^(αt). Sign of α decides spiral-in, center, or spiral-out."

Coming next: Section 23.06 handles the last linear case — repeated eigenvalues, where the eigenvector method gives only one solution and we have to build the second from a generalized eigenvector. After that, section 23.07 enters the nonlinear world, where eigenvalues of the linearization tell us the local geometry near every fixed point.