Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognise when a 2×2 system $\mathbf{x}' = A\mathbf{x}$ has a repeated eigenvalue by checking that the discriminant $\tau^2 - 4\,\det(A) = 0$ .
Distinguish the defective case (one eigenvector, improper node) from the non-defective case ( $A = \lambda I$ , star node) by inspecting $N = A - \lambda I$ .
Construct the missing second solution $\mathbf{x}_2(t) = e^{\lambda t}\,(t\,\mathbf{v} + \mathbf{w})$ using a generalised eigenvector $\mathbf{w}$ solving $(A - \lambda I)\,\mathbf{w} = \mathbf{v}$ .
Apply the closed-form matrix exponential $e^{At} = e^{\lambda t}\,(I + tN)$ when $N$ is nilpotent.
Draw the phase portrait of an improper node and explain why trajectories curl into the eigenvector line in forward time.
Cross-check by hand, by NumPy, by SciPy's general $\texttt{expm}$ , and by PyTorch's $\texttt{matrix\_exp}$ .

Why a Repeated Eigenvalue Is a Problem

In Sections 23.04 and 23.05 you saw the pleasant cases. With two distinct real eigenvalues $\lambda_1 \neq \lambda_2$ the system has two linearly independent eigenvectors $\mathbf{v}_1, \mathbf{v}_2$ , and the general solution is the textbook combination

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

With complex conjugate eigenvalues $\lambda = \alpha \pm i\beta$ , the same recipe works after Euler's formula turns the complex exponential into $e^{\alpha t}(\cos\beta t \pm i\sin\beta t)$ — a damped or growing spiral.

Now collapse the two real eigenvalues into a single number. Try

A = \begin{pmatrix} 3 & 1 \\ 0 & 3 \end{pmatrix}.

The characteristic polynomial is $(3 - r)^2 = 0$ , so the only eigenvalue is $\lambda = 3$ , repeated. Looking for eigenvectors:

(A - 3I)\,\mathbf{v} = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix} v_1 \\ v_2 \end{pmatrix} = \mathbf{0} \implies v_2 = 0.

The eigenvector line is the entire kernel. It is one- dimensional, spanned only by $\mathbf{v} = (1, 0)$ .

The trouble

The state space $\mathbb{R}^2$ is two-dimensional, but we have found only one linearly independent solution $\mathbf{x}_1(t) = e^{3t}\,\mathbf{v}$ . To match two initial conditions $x_1(0), x_2(0)$ we need two linearly independent solutions. The eigenvalue method, as written, is broken. We need a second solution out of thin air.

Intuition: A Razor's Edge Between Node and Spiral

Before we do the algebra, look at the picture. The classification of 2×2 linear systems lives in the $(\tau, \det)$ plane:

Condition on discriminant Δ = τ² − 4·det	Eigenvalues	Phase portrait
Δ > 0, det > 0, τ < 0	two distinct negative reals	Stable node — straight-line attractor
Δ > 0, det < 0	one positive, one negative real	Saddle — hyperbolic flow
Δ < 0, τ < 0	complex conjugates with negative real part	Stable spiral
Δ = 0	single repeated real eigenvalue	Improper node (defective) OR star node (non-defective) — this section
Δ > 0, det > 0, τ > 0	two distinct positive reals	Unstable node
Δ < 0, τ > 0	complex conjugates with positive real part	Unstable spiral

The repeated-eigenvalue case lives exactly on the boundary between “two real eigenvalues” and “two complex eigenvalues”. It is a measure-zero slice of parameter space — a razor's edge. Move $\Delta$ a hair into the positive half and the picture is a node with two distinct asymptotic directions. Move $\Delta$ a hair into the negative half and you get a spiral that turns through every angle. The repeated- eigenvalue picture has to do both jobs at once: it must look like a node (one direction is asymptotically preferred) and a proto-spiral (every other trajectory has to bend its way into that direction). That bending is exactly the role of the linear-in- $t$ term we are about to derive.

A picture-only summary

A real node has two preferred lines. A spiral has none. The improper node has one, and every trajectory eventually flattens onto it. The trajectory does not enter the line tangentially the way a star node does — it slides in parallel, the way a car drifts up to the curb when you turn the wheel.

Where the Missing Second Solution Comes From

Guess a second solution of the form

\mathbf{x}_2(t) = e^{\lambda t}\,\bigl(t\,\mathbf{v} + \mathbf{w}\bigr),

where $\mathbf{w}$ is some unknown constant vector and $\mathbf{v}$ is the eigenvector we already have. Why this guess? Because $t\,e^{\lambda t}$ is the second linearly independent solution of the scalar equation $(D - \lambda)^2 y = 0$ from Chapter 22 — and that is exactly the structure we are now lifting to two dimensions.

Differentiate the guess:

\mathbf{x}_2'(t) = \lambda\,e^{\lambda t}\,(t\,\mathbf{v} + \mathbf{w}) + e^{\lambda t}\,\mathbf{v}.

Plug it into $\mathbf{x}_2' = A\,\mathbf{x}_2$ :

\lambda\,(t\,\mathbf{v} + \mathbf{w}) + \mathbf{v} = A\,(t\,\mathbf{v} + \mathbf{w}).

Match coefficients of $t$ on each side. The $t$ -coefficient gives $\lambda\,\mathbf{v} = A\,\mathbf{v}$ , which is just the eigenvalue equation we already used. The constant term gives the new requirement:

\boxed{\,(A - \lambda I)\,\mathbf{w} = \mathbf{v}\,.}

That is one linear system in $\mathbf{w}$ . It has a solution because $\mathbf{v}$ is in the image of $A - \lambda I$ — a fact we can prove algebraically and which we'll feel geometrically in the next subsection.

Why the guess works

The forcing term $e^{\lambda t}\,\mathbf{v}$ is exactly the first solution, sitting on the eigenvector line. The system is asked “follow yourself, but shifted by $\mathbf{v}$ .” The only way to satisfy that without falling onto the eigenvector line is to ride a linearly-growing copy of $\mathbf{v}$ attached to an offset vector $\mathbf{w}$ . The extra factor of $t$ is the system's response to resonance with its own eigenmode.

The Generalised Eigenvector (the Algebraic Story)

A vector $\mathbf{w}$ satisfying $(A - \lambda I)\,\mathbf{w} = \mathbf{v}$ is called a generalised eigenvector of rank 2. Apply $(A - \lambda I)$ a second time and you get

(A - \lambda I)^2\,\mathbf{w} = (A - \lambda I)\,\mathbf{v} = \mathbf{0}.

So $\mathbf{w}$ lives in the kernel of $(A - \lambda I)^2$ but not in the kernel of $A - \lambda I$ itself. The chain looks like this:

\mathbf{w} \xrightarrow{A - \lambda I} \mathbf{v} \xrightarrow{A - \lambda I} \mathbf{0}.

For 2×2 matrices this two-step ladder is the entire Jordan decomposition. The pair $\{\mathbf{v}, \mathbf{w}\}$ is a basis of $\mathbb{R}^2$ , and in that basis the matrix takes its Jordan canonical form:

J = \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix},

a so-called Jordan block. The off-diagonal $1$ is the algebraic fingerprint of the defective case — it is the entry that creates the $t$ term in the solution.

Why w always exists

We need to know that the equation $N\,\mathbf{w} = \mathbf{v}$ , where $N = A - \lambda I$ , has a solution. The algebraic argument is short: $N$ is nilpotent (we'll prove this in a moment), so its image equals its kernel for a 2×2 nilpotent of index 2. The eigenvector $\mathbf{v}$ sits in the kernel, hence in the image — meaning some $\mathbf{w}$ maps to it under $N$ .

The nilpotency of N — Cayley–Hamilton in baby form

A 2×2 matrix $A$ with repeated eigenvalue $\lambda$ has characteristic polynomial $(r - \lambda)^2$ . The Cayley–Hamilton theorem says every matrix satisfies its own characteristic polynomial:

(A - \lambda I)^2 = 0.

That is the algebraic punchline of the section. Once you know $N^2 = 0$ , the power series of $e^{Nt}$ telescopes after two terms — and that is where the $(I + tN)$ factor in the matrix exponential comes from.

Two Cases: Defective vs. Non-Defective

Having a repeated eigenvalue does not automatically mean you are stuck with only one eigenvector. The subtlety is the difference between algebraic multiplicity (the multiplicity of $\lambda$ as a root of the characteristic polynomial) and geometric multiplicity (the dimension of $\ker(A - \lambda I)$ ).

	Defective (improper node)	Non-defective (star / proper node)
Matrix shape	Genuine Jordan block, e.g. [[3,1],[0,3]]	Scalar matrix λI, e.g. [[3,0],[0,3]]
N = A − λI	Nonzero, rank 1, N² = 0	Identically zero
Algebraic mult.	2	2
Geometric mult.	1	2
Independent eigenvectors	Just one, v	Every nonzero vector
General solution	e^{λt}(c₁ v + c₂(t v + w))	e^{λt}(c₁ v₁ + c₂ v₂) for any basis
Phase portrait	Curves bending into the v-line	Straight rays from origin (star)
Common name	Improper / degenerate node	Star node / proper node

Practically: compute $N = A - \lambda I$ . If $N = 0$ , you are in the star case and every direction is preserved. Otherwise $N$ has rank 1, you are defective, and you must find a generalised eigenvector.

Almost never non-defective in practice

Real matrices coming from physics are almost surely defective when they have a repeated eigenvalue, because $A = \lambda I$ is a single point in the space of all matrices. Floating-point noise from finite-difference discretisations or measurement is enough to perturb a star node into a slightly-defective configuration — or, equally often, to split the repeated eigenvalue into a near-pair.

The Matrix Exponential View

The cleanest, most uniform way to encode every case at once is the matrix exponential:

\mathbf{x}(t) = e^{A t}\,\mathbf{x}(0), \quad e^{At} = \sum_{k=0}^{\infty} \frac{(At)^k}{k!}.

Decompose $A = \lambda I + N$ . Because $\lambda I$ commutes with everything,

e^{At} = e^{(\lambda I)t}\,e^{Nt} = e^{\lambda t}\,e^{Nt}.

Now expand $e^{Nt}$ as a power series. Because $N^2 = 0$ , every term with $k \ge 2$ vanishes:

e^{Nt} = I + tN + \tfrac{1}{2!}(tN)^2 + \cdots = I + tN.

\boxed{\;e^{At} = e^{\lambda t}\bigl(I + tN\bigr)\;}

and the trajectory is

\mathbf{x}(t) = e^{\lambda t}\bigl(I + tN\bigr)\,\mathbf{x}(0).

Multiplying out gives back the $(c_1 + c_2 t)\mathbf{v} + c_2\mathbf{w}$ form from before — the two derivations are the same identity in different clothes.

The series picture

For a generic $A$ , the series $e^{At}$ is infinite. For a 2×2 matrix with a repeated eigenvalue it is just two terms. That brevity is what makes the defective case so algebraically clean once you stop fighting it.

What the Phase Portrait Actually Looks Like

Pick a generic initial condition $\mathbf{x}(0)$ that is not on the eigenvector line. Decompose it as $\mathbf{x}(0) = c_1\,\mathbf{v} + c_2\,\mathbf{w}$ with $c_2 \neq 0$ . Then

\mathbf{x}(t) = e^{\lambda t}\bigl[(c_1 + c_2 t)\,\mathbf{v} + c_2\,\mathbf{w}\bigr].

Look at the bracket. The component along $\mathbf{v}$ grows linearly in $t$ (because of the $c_2 t$ ), while the component along $\mathbf{w}$ stays fixed at $c_2$ . As $t \to \infty$ the v-component dominates, and the unit vector pointing in the direction of $\mathbf{x}(t)$ tends to $\pm\mathbf{v}$ . Geometrically: every trajectory asymptotically aligns with the eigenvector line.

Backward in time, the same argument with $t \to -\infty$ shows the v-component still wins (its magnitude grows like $|c_2 t|$ ), so trajectories also align with the eigenvector line at the other end. The trajectory is two opposite asymptotes glued together by a single bend.

The slogan for an improper node

“Every road eventually merges onto the eigenvector highway.” The merge happens parallel to $\mathbf{v}$ , not tangentially — because the offset along $\mathbf{w}$ never changes, it just becomes relatively smaller compared to the growing $v$ -component.

Direction of curvature is fixed by the sign of c₂

Initial conditions with $c_2 > 0$ bend in one direction; with $c_2 < 0$ they bend in the opposite direction. The eigenvector line itself is the trajectory with $c_2 = 0$ : it grows or decays straight, with no bending at all.

Worked Example (Step-by-Step)

Solve

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

Click to expand the full hand calculation

Step 1 — Trace and determinant.

$\tau = 3 + 3 = 6$ , $\det = 3\cdot 3 - 1\cdot 0 = 9$ . Discriminant $\Delta = \tau^2 - 4\det = 36 - 36 = 0$ — repeated eigenvalue confirmed.

Step 2 — Eigenvalue from the trace.

$\lambda = \tau / 2 = 3$ .

Step 3 — Form N = A − λI.

N = A - 3I = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}.

Check nilpotency: $N^2 = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}$ . Confirmed.

Step 4 — Eigenvector v from ker(N).

$N\,\mathbf{v} = \mathbf{0}$ : the first row says $v_2 = 0$ ; the second row is automatic. So $\mathbf{v} = (1, 0)$ .

Step 5 — Generalised eigenvector w from Nw = v.

$\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}\begin{pmatrix} w_1 \\ w_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ gives $w_2 = 1$ with $w_1$ free. Pick the simplest: $\mathbf{w} = (0, 1)$ .

Step 6 — Coefficients c₁, c₂ from the initial condition.

$\mathbf{x}(0) = c_1\,\mathbf{v} + c_2\,\mathbf{w} = (c_1, c_2) = (1, 2)$ , so $c_1 = 1,\; c_2 = 2$ .

Step 7 — Assemble the closed form.

\mathbf{x}(t) = e^{3t}\bigl[(1 + 2t)\,\mathbf{v} + 2\,\mathbf{w}\bigr] = e^{3t}\begin{pmatrix} 1 + 2t \\ 2 \end{pmatrix}.

Step 8 — Sanity-check at the boundary.

At $t = 0$ : $\mathbf{x}(0) = (1, 2)$ . ✓ Differentiate: $\mathbf{x}'(t) = 3e^{3t}(1 + 2t, 2) + e^{3t}(2, 0) = e^{3t}(3(1 + 2t) + 2,\; 6)$ . Compare with $A\,\mathbf{x}(t) = e^{3t}\begin{pmatrix} 3(1 + 2t) + 2 \\ 6 \end{pmatrix}$ . Identical. ✓

Step 9 — Evaluate at three times.

$\mathbf{x}(0.5) = e^{1.5}(2, 2) \approx (8.963,\; 8.963)$ . Notice both components are equal — the trajectory passes through the line $y = x$ at exactly this time, because $1 + 2t = 2 =$ second component.

$\mathbf{x}(1.0) = e^{3}(3, 2) \approx (60.257,\; 40.171)$ . The ratio of components is $3/2$ ; as $t \to \infty$ this ratio diverges and the trajectory aligns with $\mathbf{v} = (1, 0)$ (the x-axis).

Step 10 — Read the geometry.

The eigenvector line is the x-axis. The trajectory starts at $(1, 2)$ , has constant $w$ -component $c_2 = 2$ (the y-component stays equal to $2\,e^{3t}$ , scaled but not changing its sign), and the x-component grows like $(1 + 2t)\,e^{3t}$ — faster than the y. So the trajectory races off to the right, asymptotically parallel to the x-axis.

Interactive Phase-Portrait Explorer

Click anywhere on the phase plane to place a new initial condition. Drag the matrix entries below — the “snap d” checkbox keeps the discriminant locked to zero so you stay in the repeated- eigenvalue regime as you scrub. Toggle it off and watch the picture de-tune.

Loading repeated-eigenvalue phase portrait…

Interactive: Walk the Discriminant Through Zero

Here the matrix is $A(s) = \begin{pmatrix} -0.4 & 1 \\ s & -0.4 \end{pmatrix}$ and the discriminant is $4s$ . Sliding $s$ through zero takes you through three qualitatively different phase portraits in a single continuous motion.

Loading discriminant bifurcation viewer…

Computation in Python

The closed-form solver (no ODE integrator needed)

This routine takes $A, \mathbf{x}(0), t$ and returns the exact trajectory plus a diagnostics dictionary. It handles both the defective and non-defective cases without branching the user-facing API.

Closed-form solver for a 2×2 system with a repeated eigenvalue

🐍repeated_eigen_solver.py

Explanation(15)

Code(70)

1Import NumPy

We need eye, trace, exp, broadcasting, and pinv. The entire solver is one vectorised NumPy expression — no Python loop over time, no ODE integrator. Repeated-eigenvalue problems have a closed form; you should never integrate them numerically.

3Function header — three arguments

A is the 2×2 system matrix. x0 is the initial state (2-vector). t is a NumPy array of time samples we want the trajectory evaluated at. We return the trajectory plus a diagnostics dict so the caller can see lambda, v, w, and whether the matrix turned out to be defective.

14Coerce inputs to float arrays

Even if the user passes Python ints, we promote to float64 before doing any algebra. Integer arrays cause overflow on exp(λt) for large λt — float64 quietly handles values up to ~1e308.

18Eigenvalue from the trace — no QR, no characteristic polynomial

For a 2×2 with a repeated eigenvalue, the trace tells you everything: τ = a + d and both eigenvalues equal τ/2. We do NOT call np.linalg.eig — round-off in eig can split a true repeated eigenvalue into a near-pair of complex conjugates, which then sends the code down the wrong branch. Reading λ directly off the trace is exact.

22Build N = A − λI — the nilpotent piece

This is the algebraic heart of the section. If λ is a repeated eigenvalue, then N = A − λI satisfies N² = 0 (a 2×2 matrix with a repeated eigenvalue is at most nilpotent of index 2 by Cayley–Hamilton). N captures EXACTLY what is left over after you peel off the scalar part λI. The size of N tells you whether A is defective.

24Branch 1 — non-defective: N is the zero matrix

If N = 0 then A = λI: a pure scalar multiple of the identity. Every vector is an eigenvector. The flow is x(t) = e^{λt} x(0), every trajectory a straight ray from the origin. This is the 'star node' / 'proper node' case. Geometrically: total isotropy; algebraically: two independent eigenvectors for one eigenvalue (geometric multiplicity = algebraic multiplicity = 2).

26Vectorised x(t) = e^{λt} x(0)

exp(λ * t) is a 1-D array of shape (T,). We add a new axis to broadcast against x0 (shape (2,)), producing X of shape (T, 2). NumPy expands automatically; no Python loop.

30Branch 2 — defective: find one eigenvector v

We are now in the case where N has rank 1 (its kernel is 1-dimensional). Pick whichever row of N has the larger pivot for numerical stability, then the kernel of a [p q] row is spanned by (−q, p). Geometric multiplicity is 1 while algebraic multiplicity is 2 — this is the 'defective' or 'deficient' or 'improper-node' regime.

35Normalise v

Cosmetic — eigenvectors are defined up to scale. Normalising keeps the coefficients c1, c2 in a human-readable range and avoids any one component from dominating the others purely because of magnitude.

38Generalised eigenvector w via pinv

We need ANY solution w to N w = v. Such a w exists precisely because N is nilpotent: v ∈ image(N) = kernel(N). The Moore–Penrose pseudoinverse pinv(N) returns the minimum-norm solution. (We could also row-reduce by hand, but pinv is a one-liner that always works.)

41Express x0 in the (v, w) basis

v and w are linearly independent (provably — if they were collinear, w would be in kernel(N), and then N w = 0 ≠ v). So they form a basis of ℝ². np.linalg.solve hands back the unique coefficients (c1, c2) such that x0 = c1·v + c2·w. The c2 coefficient is the one that 'turns on' the t·v term.

45Closed form — the famous (c1 + c2 t) v + c2 w combination

This is the entire section in one expression. The envelope e^{λt} handles growth/decay. Inside, the term (c1 + c2 t) v rides the eigenvector direction at a linearly-growing magnitude — THIS is where the extra t comes from. The c2·w term anchors the trajectory off the eigenvector line so it has somewhere to curve back from.

49Vectorised broadcast over the entire time grid

env has shape (T,). (c1 + c2 t)[:, None] has shape (T, 1) so it broadcasts against v (shape (2,)) to give a (T, 2) array. We add the constant offset c2 * w (shape (2,)) and finally multiply by the (T, 1) envelope. The whole trajectory is one expression — no Python-level loop.

53Worked example: A = [[3,1],[0,3]], x(0) = (1, 2)

Trace = 6, so λ = 3. N = A − 3I = [[0,1],[0,0]], so kernel(N) is spanned by (1, 0) — that is v. Solving N w = v gives w = (0, 1). x(0) = (1, 2) = 1·v + 2·w, so c1 = 1, c2 = 2. The closed form is x(t) = e^{3t} (1 + 2t, 2).

60Print diagnostics

λ = 3.0 (exact from trace), v = (1, 0), w = (0, 1). At t = 1: x(1) = e³ · (3, 2) ≈ (60.257, 40.171). At t = 0.5: x(0.5) = e^{1.5} · (2, 2) ≈ (8.963, 8.963) — note BOTH components equal e^{1.5}·2 because the first component is 1 + 2·0.5 = 2, identical to the second. The trajectory hits the line y = x at exactly t = 0.5.

55 lines without explanation

1import numpy as np
2
3def solve_repeated(A, x0, t):
4    """
5    Solve x'(t) = A x  for a 2x2 matrix A that has a REPEATED real eigenvalue.
6
7    Handles both sub-cases:
8      * defective  (only one independent eigenvector v)
9      * non-defective  (A = lambda * I, every direction is an eigenvector)
10
11    Returns
12    -------
13    X : np.ndarray of shape (len(t), 2)   the trajectory.
14    info : dict with the eigenvalue, eigenvector v, generalised eigenvector w,
15           and whether A is defective.
16    """
17    A = np.asarray(A, dtype=float)
18    x0 = np.asarray(x0, dtype=float)
19
20    # 1. Repeated eigenvalue from the trace:  lambda = trace / 2
21    lam = 0.5 * np.trace(A)
22
23    # 2. Look at  N = A - lambda I.   N must be NILPOTENT  (N**2 = 0)
24    #    when lambda is repeated.  This is the algebraic source of the t-term.
25    N = A - lam * np.eye(2)
26
27    if np.allclose(N, 0.0):
28        # Non-defective:  A = lambda I.   Pure radial expansion / contraction.
29        X = np.exp(lam * t)[:, None] * x0[None, :]
30        return X, {"lambda": lam, "v": None, "w": None, "defective": False}
31
32    # 3. Defective case.   Find ONE eigenvector v in kernel of N.
33    #    Use the row with the largest pivot.
34    if abs(N[0, 0]) + abs(N[0, 1]) >= abs(N[1, 0]) + abs(N[1, 1]):
35        v = np.array([-N[0, 1], N[0, 0]])
36    else:
37        v = np.array([-N[1, 1], N[1, 0]])
38    v = v / np.linalg.norm(v)
39
40    # 4. Find generalised eigenvector w with   N w = v   (any solution).
41    #    Use pinv:  it is the minimum-norm exact solution.
42    w = np.linalg.pinv(N) @ v
43
44    # 5. Decompose x0 in the (v, w) basis:   x0 = c1 v + c2 w
45    M = np.column_stack([v, w])              # 2x2 invertible
46    c1, c2 = np.linalg.solve(M, x0)
47
48    # 6. Closed form:  x(t) = e^{lam t} [ (c1 + c2 t) v + c2 w ]
49    env = np.exp(lam * t)
50    coef = (c1 + c2 * t)[:, None] * v[None, :]   # shape (T, 2)
51    extra = c2 * w[None, :]                       # broadcast over time
52    X = env[:, None] * (coef + extra)
53    return X, {"lambda": lam, "v": v, "w": w, "defective": True}
54
55
56# ---- run the worked example ----
57A = np.array([[3.0, 1.0],
58              [0.0, 3.0]])
59x0 = np.array([1.0, 2.0])
60t = np.linspace(0.0, 1.0, 6)
61
62X, info = solve_repeated(A, x0, t)
63
64print("lambda     =", info["lambda"])      # 3.0
65print("v          =", info["v"])           # (1, 0)
66print("w          =", info["w"])           # (0, 1)
67print("defective? =", info["defective"])   # True
68print("x(0.0)     =", X[0])                # (1, 2)
69print("x(0.5)     =", X[2])                # ≈ (8.963, 8.963)
70print("x(1.0)     =", X[-1])               # ≈ (60.257, 40.171)

Verifying with the matrix exponential

The same trajectory can be computed by exponentiating $A$ directly. For a Jordan block this gives the compact closed form $e^{At} = e^{\lambda t}(I + tN)$ — and we cross-check against $\texttt{scipy.linalg.expm}$ .

Verifying e^{At} = e^{λt}(I + tN) for a 2×2 Jordan block

🐍repeated_eigen_matrix_exp.py

Explanation(8)

Code(26)

2Import scipy.linalg.expm

scipy.linalg.expm evaluates the matrix exponential e^M for any square M using Padé approximation with scaling-and-squaring. We use it ONLY as a sanity-check reference — the closed form we derived is exact for this particular A.

8Eigenvalue from the trace

Exactly as in the previous solver: τ = a + d = 6, so λ = 3. No solver call needed for a repeated eigenvalue.

9The nilpotent piece N = A − λI

Computing N is one matrix subtraction. For A = [[3,1],[0,3]] and λ = 3, this gives N = [[0,1],[0,0]] — the standard 2×2 nilpotent matrix that maps the y-axis to the x-axis and squashes everything else.

12Assert N is nilpotent

Algebraic fact: for any 2×2 matrix with a repeated eigenvalue λ, (A − λI)² = 0. This is a baby case of the Cayley–Hamilton theorem (every matrix satisfies its own characteristic polynomial, here (r − λ)² = 0). The assert documents this invariant in the code.

16Closed-form matrix exponential of a Jordan block

Because λI and N commute (λI commutes with everything) AND N² = 0, the series e^{At} = Σ (At)^k / k! telescopes after only TWO terms in the N-direction: e^{At} = e^{λt} · e^{Nt} = e^{λt} (I + tN + 0 + 0 + …) = e^{λt} (I + tN). This is THE source of the extra t-factor that appears in the general solution.

17Function expA(t)

Returns the 2×2 matrix e^{At} evaluated at one scalar time t. Multiplying by x(0) on the right gives the state at that t. No iteration, no solver — closed form.

21Cross-check against scipy

We compute e^{A·1} two different ways: our closed form and scipy's general-purpose Padé approximation. They MUST agree to ~1e-12. If they don't, our derivation is wrong. (Spoiler: they agree — both give [[e³, e³], [0, e³]] ≈ [[20.086, 20.086], [0, 20.086]].)

26Evaluate the trajectory at three sample times

expA(t) @ x0 maps the initial state through e^{At}. For t = 0 we get (1, 2) trivially. For t = 0.5 we get e^{1.5} · (1 + 0.5·2, 2) = e^{1.5} · (2, 2) ≈ (8.963, 8.963). For t = 1 we get e³ · (1 + 1·2, 2) = e³ · (3, 2) ≈ (60.257, 40.171). Same numbers as the previous solver — both routes converge.

18 lines without explanation

1import numpy as np
2from scipy.linalg import expm   # matrix exponential
3
4A = np.array([[3.0, 1.0],
5              [0.0, 3.0]])
6
7# Decompose A = lambda I + N   (Jordan form for repeated eigenvalues)
8lam = 0.5 * np.trace(A)        # 3.0
9N   = A - lam * np.eye(2)      # [[0, 1], [0, 0]]
10
11# Sanity check the nilpotency:  N @ N must be the zero matrix.
12assert np.allclose(N @ N, 0.0), "N should be nilpotent"
13
14# Closed-form matrix exponential for a Jordan block:
15#     e^{At} = e^{lambda t} (I + t N)
16def expA(t):
17    return np.exp(lam * t) * (np.eye(2) + t * N)
18
19# Cross-check at t = 1 against scipy's general-purpose matrix exponential.
20print("closed form e^{A·1} =\n", expA(1.0))
21print("scipy expm(A·1)     =\n", expm(A * 1.0))
22
23# Trajectory from x0 = (1, 2)
24x0 = np.array([1.0, 2.0])
25for t in [0.0, 0.5, 1.0]:
26    print(f"x({t}) = {expA(t) @ x0}")

PyTorch View: matrix_exp Handles Every Case

Once you commit to the matrix-exponential picture, there is no longer a special-case branch for repeated eigenvalues. The same line of code $X = e^{At} X(0)$ works for every regime, on GPUs, batched, and differentiable through the entries of $A$ . That is the modern way the recurring building block of linear control theory enters deep learning (Linear Neural ODEs, S4, Mamba state-space models).

Batched, differentiable trajectory via torch.linalg.matrix_exp

🐍repeated_eigen_torch.py

Explanation(8)

Code(27)

1Import torch

We need torch.linalg.matrix_exp (added in 1.7) plus broadcasting and einsum. The big payoff over NumPy is that A can be made trainable: every operation below is differentiable, so gradient-based learning of dynamical systems falls out for free.

4A as a tensor, x0 as a tensor

Identical matrix and initial state as the NumPy run. requires_grad=False here because we only want a forward pass. Flip it to True and torch will accumulate ∂x(t)/∂A through the entire chain that follows.

9Time grid as a tensor

T is a 1-D tensor of 6 evenly-spaced sample times in [0, 1]. We will evaluate the WHOLE trajectory in parallel — one matrix exponential per time, vectorised on GPU if available.

12Broadcast t against A to get a batch of A*t matrices

T[:, None, None] has shape (6, 1, 1). Multiplying by A (shape (2, 2)) broadcasts up to (6, 2, 2). Each slice At[k] is the matrix A · T[k]. No Python loop.

16Batched matrix exponential

torch.linalg.matrix_exp computes e^M for every 2×2 slice in parallel. Crucially, it handles every regime — distinct real eigenvalues, complex conjugate eigenvalues, AND defective repeated eigenvalues — without any special-case code on our side. The result is a (6, 2, 2) tensor where slice k is e^{A · T[k]}.

20einsum 'kij,j->ki' — batch matrix-vector product

For each k, multiply the 2×2 matrix expAt[k] by the 2-vector x0. The output X[k] is the state at time T[k]. einsum lets us write the index pattern directly — much clearer than expAt @ x0[None, :] gymnastics. X has shape (6, 2).

23Print one row per time

Same numbers as the NumPy and SciPy paths: x(0) = (1, 2), x(0.2) ≈ e^{0.6} · (1.4, 2) ≈ (2.552, 3.644), …, x(1.0) ≈ (60.257, 40.171). All three routes agree to machine precision.

27Why PyTorch and not just SciPy?

Three reasons. First, autodiff: with A.requires_grad_(True) you can backprop through the entire trajectory and FIT a matrix A to observed data — exactly what neural-ODE and structured-state-space models (S4, Mamba) do internally. Second, GPU batching: simulating thousands of (A, x0) pairs in parallel for parameter sweeps. Third, this is the building block of linear control-theory layers — the eigenstructure of A IS the inductive bias of the model.

19 lines without explanation

1import torch
2
3# Same matrix A and initial condition as before, but as PyTorch tensors.
4A  = torch.tensor([[3.0, 1.0],
5                   [0.0, 3.0]], requires_grad=False)
6x0 = torch.tensor([1.0, 2.0])
7
8# Evaluate trajectory at many times in PARALLEL.
9T  = torch.linspace(0.0, 1.0, 6)             # shape (6,)
10
11# Build a batch of (T, 2, 2) matrices A*t, one per time sample.
12At = T[:, None, None] * A                    # broadcasting: (6, 1, 1) * (2, 2)
13
14# torch.linalg.matrix_exp handles every Jordan-block case correctly,
15# including the repeated-eigenvalue defective case we are studying.
16expAt = torch.linalg.matrix_exp(At)          # shape (6, 2, 2)
17
18# Apply each e^{A t_k} to x0.  einsum is the cleanest expression:
19#    X[k, i] = sum_j  expAt[k, i, j] * x0[j]
20X = torch.einsum("kij,j->ki", expAt, x0)     # shape (6, 2)
21
22for k, t in enumerate(T):
23    print(f"x({t.item():.2f}) = {X[k].tolist()}")
24
25# Why bother with PyTorch?  Because A can have requires_grad=True and we
26# can differentiate the trajectory through A's entries — the basis of every
27# neural-ODE / state-space model training loop.

Why this matters for ML

A trainable matrix $A$ can pass through a defective configuration during training, even if it started generic. The matrix exponential view doesn't care — gradients flow continuously through that razor's edge. Algorithms that branch on “is $A$ diagonalisable?” suddenly produce non-differentiable training surfaces and crash.

Where Repeated Eigenvalues Actually Show Up

⚙️ Critically damped mechanical systems

The mass-spring-damper with $c^2 = 4mk$ sits at the boundary between oscillating and over-damped. The state-space matrix has a single repeated eigenvalue $\lambda = -c/(2m)$ .
Designers aim for this regime because it is the fastest non-oscillating return to rest — automotive shock absorbers, door-closers, gun recoil mechanisms.

⚡ Critically damped circuits

Series RLC with $R^2 = 4L/C$ has a repeated eigenvalue — the fastest settling of a power-supply transient.
Audio-amplifier compensation networks, oscilloscope-probe compensation, snubber circuits.

🎯 Control theory

Pole placement: engineers deliberately assign repeated closed-loop poles to a plant for minimum-time non- overshooting step response.
State observers (Luenberger) often use repeated eigenvalues for the observer poles to keep the error dynamics tight.

🧠 Machine learning

S4 and Mamba state-space layers parameterise $e^{At}$ directly. Repeated eigenvalues are a measure-zero event in the parameter manifold — they are the points where eigendecomposition is discontinuous, so models that work in eigen-basis break.
The matrix-exponential parameterisation is precisely what makes those models smooth across the defective ridge.

Common Pitfalls

Don't forget the t term

A frequent student mistake in the defective case is to write the general solution as $c_1\,e^{\lambda t}\,\mathbf{v} + c_2\,e^{\lambda t}\,\mathbf{w}$ without the $t\,\mathbf{v}$ piece. That expression is not a solution of $\mathbf{x}' = A\mathbf{x}$ — you can check by differentiation that the residual is $-c_2\,e^{\lambda t}\,\mathbf{v}$ , which is nonzero unless $c_2 = 0$ .

Generalised eigenvector is not unique

Any $\mathbf{w} + s\,\mathbf{v}$ with $s \in \mathbb{R}$ is also a valid generalised eigenvector, because $N(\mathbf{w} + s\mathbf{v}) = \mathbf{v} + 0 = \mathbf{v}$ . Different choices simply re-parameterise $c_1$ . Pick whatever $\mathbf{w}$ is cheapest to compute.

np.linalg.eig may split the repeated eigenvalue

Floating-point eigensolvers detect the defective case poorly. Calling $\texttt{np.linalg.eig}([[3, 1], [10^{-12}, 3]])$ returns two distinct “eigenvalues” differing by $\sim 10^{-6}$ and two “eigenvectors” that point in nearly the same direction. Trust the trace, not the solver, when you know $\Delta = 0$ .

Star node is not improper

If you ever land on $A = \lambda I$ , do NOT compute a generalised eigenvector. There is no defective structure to repair: every direction is already an eigenvector, and the trajectories are straight rays. Confusing the two is the most common conceptual error in the section.

Summary

A 2×2 linear system $\mathbf{x}' = A\mathbf{x}$ with discriminant $\Delta = \tau^2 - 4\det A = 0$ has a single repeated eigenvalue $\lambda = \tau/2$ . The recipe is:

Compute $N = A - \lambda I$ . If $N = 0$ you are in the non-defective star case — solution is $e^{\lambda t}\,\mathbf{x}(0)$ with no further work.
Otherwise find $\mathbf{v} \in \ker N$ .
Solve $N\,\mathbf{w} = \mathbf{v}$ for a generalised eigenvector $\mathbf{w}$ .
Decompose $\mathbf{x}(0) = c_1\,\mathbf{v} + c_2\,\mathbf{w}$ .
Write $\mathbf{x}(t) = e^{\lambda t}\bigl[(c_1 + c_2 t)\,\mathbf{v} + c_2\,\mathbf{w}\bigr]$ , equivalently $e^{\lambda t}(I + tN)\,\mathbf{x}(0)$ .

Where it sits in the calculus of dynamical systems

Δ = τ² − 4·det	Eigenvalues	Section	Picture
> 0	distinct real	23.04	Node / saddle
< 0	complex conjugate	23.05	Spiral / centre
= 0 (this section)	single repeated real	23.06	Improper or star node
any	general n × n	23.07 →	Linearisation around fixed points of nonlinear systems

The slogan

“When the eigenvectors run out, the matrix exponential steps in.

e^{At} = e^{\lambda t}(I + tN)

is the entire chapter in one expression.”

Coming next: Real linear systems are an idealisation. Most physical, biological, and economic models are nonlinear. In Section 23.07 we will see that the eigenvalue machinery you just mastered — including the repeated-eigenvalue razor's edge — survives almost intact, because near any equilibrium a nonlinear system looks like a linear one, and the classification of equilibria is decided by the same $(\tau, \det)$ plane.