Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Detect the repeated-root case by the discriminant $\Delta = b^2 - 4ac = 0$ and identify the unique root $r = -b/(2a)$ .
Explain three independent reasons why $t\,e^{rt}$ is the missing second solution: a limit of two close exponentials, reduction of order, and direct verification.
Write the general solution $y(t) = (C_1 + C_2 t)\,e^{rt}$ and apply two initial conditions to solve for $C_1, C_2$ .
Recognise this case as the boundary between overdamped and underdamped behaviour — the regime engineers call critical damping.
Implement an analytic solver, an RK4 numerical solver, and a PyTorch matrix-exponential solver, and verify that all three agree to machine precision.
See the deeper linear-algebra reason: the companion matrix has a Jordan block, and the matrix exponential of a Jordan block naturally carries a factor of $t$ .

The Puzzle: Only One Exponential

We learned in Section 22.01 that every homogeneous second-order linear ODE with constant coefficients

a y'' + b y' + c y = 0,\quad a \neq 0

is solved by trying $y = e^{rt}$ , which reduces the differential equation to the algebraic characteristic equation $a r^2 + b r + c = 0$ . The quadratic has two roots, and those roots become the two independent exponential modes of the system.

But there is a sharp edge case. What happens when $\Delta = b^2 - 4ac$ is exactly $0$ ?

r = \frac{-b \pm \sqrt{0}}{2a} = -\frac{b}{2a} \quad \text{(once)}

We get only one exponent. The naive recipe would hand us a one-parameter family of solutions $y = C\,e^{rt}$ . But a second-order ODE is supposed to need two initial conditions and have a two-parameter solution family. Something is missing.

The headline

When the characteristic equation has a single repeated root $r = -b/(2a)$ , the second linearly independent solution is $y_2(t) = t\,e^{rt}$ , and the general solution is

\boxed{\,y(t) = (C_1 + C_2\,t)\,e^{rt}\,}

The rest of this section explains, from three independent angles, where that mysterious factor of $t$ actually comes from.

Where Does t Come From? The Coalescence Argument

The cleanest derivation is to perturb the problem. Imagine the characteristic equation does not quite have a repeated root — instead the two roots are $r_1 = r - \varepsilon$ and $r_2 = r + \varepsilon$ for some tiny $\varepsilon > 0$ . Then we are in the distinct-real-roots case, where two known independent solutions are

y_1(t) = e^{(r-\varepsilon)t},\qquad y_2(t) = e^{(r+\varepsilon)t}.

Because the ODE is linear, any linear combination of these is still a solution. Consider the particular combination

Y_\varepsilon(t) \;=\; \frac{e^{(r+\varepsilon)t} - e^{(r-\varepsilon)t}}{2\varepsilon}.

This is just a scaled divided difference — the symmetric Newton quotient of $e^{st}$ with respect to the parameter $s$ , sampled at $s = r$ . Now let $\varepsilon \to 0$ .

Factor out $e^{rt}$ :

Y_\varepsilon(t) = e^{rt}\,\frac{e^{\varepsilon t} - e^{-\varepsilon t}}{2\varepsilon} = e^{rt}\,\frac{\sinh(\varepsilon t)}{\varepsilon}.

Use the Taylor expansion $\sinh(\varepsilon t) = \varepsilon t + \tfrac{(\varepsilon t)^3}{6} + \cdots$ :

\frac{\sinh(\varepsilon t)}{\varepsilon} = t + \tfrac{\varepsilon^2 t^3}{6} + \cdots \;\xrightarrow{\varepsilon \to 0}\; t.

Therefore

\lim_{\varepsilon \to 0} Y_\varepsilon(t) \;=\; t\,e^{rt}.

What this argument is really saying

As the two roots collide on the real axis, one of the two independent solutions does not vanish — it rotates into a derivative-with-respect-to-the-root. Concretely,

t\,e^{rt} \;=\; \left.\frac{\partial}{\partial s} e^{st}\right|_{s=r}.

The partial derivative with respect to the eigenvalue is what plays the role of a second independent solution when the eigenvalue is repeated. You will see this exact pattern again for eigenvalues of matrices, for Wronskian arguments, and for the matrix exponential later in this section.

Why is this still a solution after the limit?

Each $Y_\varepsilon$ solves a slightly different ODE (with discriminant $4\varepsilon^2 > 0$ ). But as $\varepsilon \to 0$ , the coefficients of that ODE converge to $a, b, c$ , and $Y_\varepsilon$ converges to $t e^{rt}$ uniformly on compact $t$ -intervals. By continuous dependence of linear ODE solutions on their coefficients, the limit solves the limiting equation — i.e. the repeated-root ODE. You can also verify directly by substitution, which we do in a moment.

Second Derivation: Reduction of Order

Here is a second, completely independent route to the same answer — and one that generalises far beyond constant coefficients.

Suppose we already know one solution $y_1 = e^{rt}$ . Look for a second solution of the form

y_2(t) = u(t)\,e^{rt}

where $u(t)$ is an unknown function, not just a constant. (If $u$ were constant, this would just be a multiple of $y_1$ , which is not linearly independent.)

Compute derivatives using the product rule:

y_2 = u\,e^{rt},\qquad y_2' = (u' + r u)\,e^{rt},\qquad y_2'' = (u'' + 2r u' + r^2 u)\,e^{rt}.

Substitute into $a y'' + b y' + c y = 0$ and divide by $e^{rt}$ (always positive):

a(u'' + 2r u' + r^2 u) + b(u' + r u) + c u = 0.

Group by derivatives of $u$ :

a\,u'' + (2 a r + b)\,u' + (a r^2 + b r + c)\,u = 0.

Two things now collapse beautifully:

The bracket on $u$ is $a r^2 + b r + c = 0$ — the characteristic equation, which holds because $r$ is its root.
The bracket on $u'$ is $2 a r + b$ . Plug $r = -b/(2a)$ : $2 a \cdot (-b/(2a)) + b = -b + b = 0$ .

Both brackets vanish! What remains is

a\,u'' = 0 \quad\Longrightarrow\quad u'' = 0 \quad\Longrightarrow\quad u(t) = C_1 + C_2\,t.

Therefore the most general solution is

y(t) = u(t)\,e^{rt} = (C_1 + C_2\,t)\,e^{rt}.

What just happened — the deep reason

The cancellation of the $u'$ bracket is not a coincidence. The quantity $2 a r + b$ is the derivative of the characteristic polynomial evaluated at $r$ : $p'(r) = 2 a r + b$ . A root is repeated precisely when $p(r) = 0$ AND $p'(r) = 0$ . So both brackets vanish together exactly when the root is repeated — which is exactly when we need a polynomial factor to fill the gap.

Third Derivation: Just Substitute and Check

The most elementary route is brute force. Take $y_2 = t\,e^{rt}$ as a candidate and plug in:

y_2' = (1 + r t)\,e^{rt},\qquad y_2'' = (2r + r^2 t)\,e^{rt}.

Then

a y_2'' + b y_2' + c y_2 = \bigl[\,(a r^2 + b r + c)\,t + (2 a r + b)\,\bigr] e^{rt}.

The coefficient of $t$ is zero because $r$ is a root of the characteristic equation. The constant piece is zero because $r = -b/(2a)$ forces $2ar + b = 0$ . So $y_2 = t e^{rt}$ is genuinely a solution.

And it is linearly independent from $y_1 = e^{rt}$ : the Wronskian is

W(t) = \begin{vmatrix} e^{rt} & t e^{rt} \\ r e^{rt} & (1 + r t) e^{rt} \end{vmatrix} = e^{2 r t}\bigl[(1 + r t) - r t\bigr] = e^{2 r t} \neq 0.

Since $W(t) \neq 0$ for every $t$ , the two solutions form a fundamental set and span the entire solution space.

The General Solution and Two Initial Conditions

Putting all three derivations together, the general solution of $a y'' + b y' + c y = 0$ with $\Delta = 0$ is

y(t) = (C_1 + C_2\,t)\,e^{rt},\qquad r = -\frac{b}{2a}.

To pin down the two constants we use both initial conditions. Setting $t = 0$ :

y(0) = C_1.

Differentiate $y = (C_1 + C_2 t) e^{rt}$ with the product rule:

y'(t) = C_2 e^{rt} + r(C_1 + C_2 t) e^{rt} = \bigl(C_2 + r C_1 + r C_2 t\bigr) e^{rt}.

Evaluate at $t = 0$ :

y'(0) = C_2 + r\,C_1 \quad\Longrightarrow\quad C_2 = y'(0) - r\,y(0).

Quantity	Formula	Comes from
Repeated root r	r = −b / (2a)	Discriminant is 0
C₁	C₁ = y(0)	Set t = 0 in y(t)
C₂	C₂ = y'(0) − r·y(0)	Set t = 0 in y'(t)
Number of zeros of y	0 or 1 on (0, ∞)	Linear · positive exponential

A useful qualitative fact

Because $e^{rt} > 0$ for every real $t$ , the sign of $y(t)$ is determined entirely by the linear factor $C_1 + C_2 t$ . A repeated-root solution therefore changes sign at most once, at $t = -C_1 / C_2$ (provided $C_2 \neq 0$ ). It is the most boring oscillation in the world — no oscillation at all.

Interactive: Watch the Limit Happen

The picture below has two coupled plots. In the top plot, drag the ε slider toward zero and watch the green divided-difference curve glue itself onto the amber $t\,e^{rt}$ curve. That is the second solution being born as the two roots collide.

In the bottom plot you can dial in any $r, C_1, C_2$ and see how the linear factor $C_1 + C_2 t$ rides on top of the exponential envelope. Pay attention to the red dashed line marking the zero crossing at $t = -C_1/C_2$ .

Loading interactive coalescence demo…

Three things to try

Set $r = 0$ . The envelope becomes 1 and $y(t) = C_1 + C_2 t$ is a literal straight line. The repeated-root machinery is consistent with this edge case.
Set $r < 0$ and $C_1 \cdot C_2 < 0$ . The curve overshoots zero once, then asymptotes back to zero from the other side. This is the critically damped overshoot — your car suspension hitting a pothole.
Set $r > 0$ . The polynomial factor is irrelevant in the long run; the growing exponential dominates. This is what happens when a controller goes unstable.

Physical Meaning: Critical Damping

Translate to the canonical mass–spring–damper $m\ddot x + c \dot x + k x = 0$ . The discriminant condition $\Delta = 0$ becomes

c^2 = 4 m k \;\Longleftrightarrow\; \zeta := \frac{c}{2\sqrt{mk}} = 1.

$\zeta$ is the damping ratio. The repeated-root case sits exactly at $\zeta = 1$ — the knife-edge between two very different behaviours.

Regime	ζ	Roots	What it looks like
Overdamped	ζ > 1	Two distinct negative real	Slow, non-oscillating return; energy bleeds away through two timescales.
Critically damped	ζ = 1	One repeated negative real	Fastest non-oscillating return. At most one zero crossing.
Underdamped	0 < ζ < 1	Complex conjugate pair	Decaying oscillation; system rings before settling.

Why engineers care: car suspensions, robotic-arm controllers, hard-disk read heads, and analog filter circuits are often tuned just barely underdamped of critical. Critical damping returns to equilibrium fastest without overshoot, but the slightest disturbance pushes a real system either over (sluggish) or under (ringing) — so designers leave a safety margin.

The boundary is mathematically thin, physically wide

In any continuous family of ODEs, the set $\{\Delta = 0\}$ is a one-dimensional surface in the three-dimensional $(a, b, c)$ parameter space — a thin slice. But it is the slice that separates the oscillating and non-oscillating worlds. Almost every interesting control-engineering decision lives within a small neighbourhood of this slice.

Worked Example (Step-by-Step)

Solve $y'' + 4 y' + 4 y = 0$ subject to $y(0) = 1$ and $y'(0) = -3$ .

Click to expand the full hand calculation

Step 1 — Read off a, b, c.

Comparing with $a y'' + b y' + c y = 0$ : $a = 1,\; b = 4,\; c = 4$ .

Step 2 — Discriminant.

$\Delta = b^2 - 4ac = 16 - 16 = 0$ . Repeated-root case confirmed.

Step 3 — The repeated root.

$r = -\dfrac{b}{2a} = -\dfrac{4}{2} = -2$ . Negative root ⇒ both basic solutions $e^{-2t}$ and $t\,e^{-2t}$ decay to zero.

Step 4 — General solution.

y(t) = (C_1 + C_2 t)\,e^{-2 t}.

Step 5 — Apply $y(0) = 1$ .

$y(0) = (C_1 + 0)\,e^{0} = C_1$ , so $C_1 = 1$ .

Step 6 — Differentiate.

y'(t) = C_2 e^{-2t} - 2(C_1 + C_2 t) e^{-2t} = \bigl(C_2 - 2 C_1 - 2 C_2 t\bigr) e^{-2t}.

Step 7 — Apply $y'(0) = -3$ .

$y'(0) = C_2 - 2 C_1 = C_2 - 2 = -3$ , so $C_2 = -1$ .

Step 8 — Particular solution.

\boxed{\,y(t) = (1 - t)\,e^{-2t}\,}

Step 9 — Sanity check the ICs.

$y(0) = 1 \cdot 1 = 1$ ✓. For $y'(0)$ : at $t = 0$ the formula $(C_2 - 2 C_1 - 2 C_2 t)e^{-2t}$ is $(-1) - 2(1) - 0 = -3$ ✓.

Step 10 — Sanity check the ODE.

Compute $y'(t) = (-1)\,e^{-2t} - 2(1 - t)\,e^{-2t} = (-3 + 2t)\,e^{-2t}$ and $y''(t) = 2\,e^{-2t} - 2(-3 + 2t)\,e^{-2t} = (8 - 4t)\,e^{-2t}$ . Then $y'' + 4 y' + 4 y$ equals $\bigl[(8 - 4t) + 4(-3 + 2t) + 4(1 - t)\bigr] e^{-2t} = (8 - 4t - 12 + 8 t + 4 - 4 t) e^{-2t} = 0$ ✓.

Step 11 — Where does y(t) cross zero?

Set $1 - t = 0 \Rightarrow t = 1$ . Because the exponential is never zero, this is the unique zero for $t > 0$ . Numerically: $y(1) = 0 \cdot e^{-2} = 0$ .

Step 12 — A long-term value.

At $t = 2$ : $y(2) = (1 - 2)\,e^{-4} = -e^{-4} \approx -0.01832$ . Check the visualization above with $r = -2,\,C_1 = 1,\,C_2 = -1$ — the curve dips below zero around $t = 1$ and returns to zero from below.

Sweep Through the Three Cases

The interactive below is the same explorer from Section 22.01, but now look at it through the lens of this section. Drag the $b$ slider slowly while watching the two root markers in the complex plane on the right. The instant they collide on the real axis is exactly the repeated-root case. Push a hair further and they pop off the real axis as a complex pair.

Loading characteristic-equation explorer…

Two presets to compare

Click Critically damped (repeated root). Note the two root markers collapse to one point. Coefficients $a = 1, b = 4, c = 4$ give $r = -2$ — exactly our worked example (the ICs differ; the algebra is the same).
Now click Overdamped (distinct real roots). The two markers split apart along the real axis. The shape of the solution is qualitatively similar (no oscillation, decays to zero), but the closed form is $C_1 e^{r_1 t} + C_2 e^{r_2 t}$ — no $t$ factor. The repeated-root case is the limit of this as $r_2 \to r_1$ .

Computation in Python — Analytic Solver

Mathematics first, then code. With the closed form $y(t) = (C_1 + C_2 t)\,e^{rt}$ in hand, a complete Python solver is short. The only judgement call is what to do when the user hands us non-repeated coefficients — we raise an explicit error rather than silently producing nonsense.

Direct evaluation of (C₁ + C₂t)·e^(rt) on a time grid

🐍repeated_root_analytic.py

Explanation(9)

Code(44)

1Import NumPy

We will evaluate the closed-form solution y(t) on a whole time grid in one shot. NumPy's exp is vectorized — multiplying a scalar by an array and exponentiating an array stays elementwise. No Python loop is needed.

3Function header

Three coefficients (a, b, c), two initial conditions (y0, yp0), and a NumPy array t of time samples. The two ICs are not optional — a second-order ODE has a 2-parameter family of solutions, so we need two pieces of starting data to pin one down.

21Discriminant guard

Δ = b² − 4ac. The repeated-root case is defined by Δ = 0 EXACTLY. We allow a tiny floating-point tolerance (1e-9) and raise otherwise — the formula below would not give two independent solutions for distinct or complex roots.

25The repeated root r = −b / (2a)

The quadratic ar² + br + c = 0 collapses to a(r − r₀)² = 0 when Δ = 0, with r₀ = −b/(2a). That single value is the only exponent that ever appears in the solution — both linearly independent solutions e^(rt) and t·e^(rt) share it.

30C₁ from y(0)

At t = 0, (C₁ + C₂·0)·e^(0) = C₁. So C₁ = y(0) directly. No system to solve yet — the second initial condition has not been used.

31C₂ from y'(0)

Differentiate y = (C₁ + C₂t)·e^(rt) using the product rule: y' = C₂·e^(rt) + r(C₁ + C₂t)·e^(rt). At t = 0 this is C₂ + r·C₁. Setting equal to y'(0) and solving: C₂ = y'(0) − r·y(0). Two ICs ⇒ two constants ⇒ unique solution.

33Build y(t) on the grid

NumPy evaluates (C₁ + C₂·t) elementwise (linear in t), then multiplies by exp(r·t) elementwise. The result is a 1-D array the same length as t. This single line IS the analytic solution.

38Run the worked example

We pick a = 1, b = 4, c = 4 — the canonical critically damped case. Δ = 16 − 16 = 0, so r = −2. With y(0) = 1 and y'(0) = −3, we expect C₁ = 1 and C₂ = −3 − (−2)(1) = −1, so y(t) = (1 − t)e^(−2t).

43Sanity print: zero crossing

Setting C₁ + C₂·t = 0 gives t = −C₁/C₂ = 1.0. That is the ONLY zero of y(t) on (0, ∞), because the exponential factor is never zero. A repeated-root solution can cross zero at most once — a sharp qualitative fact you should always check.

35 lines without explanation

1import numpy as np
2
3def solve_repeated_root(a, b, c, y0, yp0, t):
4    """
5    Solve a*y'' + b*y' + c*y = 0 for the REPEATED-ROOT case (Δ = 0).
6
7    Parameters
8    ----------
9    a, b, c : float
10        Constant coefficients. Must satisfy b**2 - 4*a*c == 0.
11    y0      : float   y(0),  initial value.
12    yp0     : float   y'(0), initial slope.
13    t       : np.ndarray, time grid on which to evaluate y.
14
15    Returns
16    -------
17    y    : np.ndarray, the same shape as t.
18    info : dict with the root r and constants C1, C2.
19    """
20    disc = b**2 - 4 * a * c
21    if abs(disc) > 1e-9:
22        raise ValueError(f"Not a repeated-root case: Δ = {disc:.4f}")
23
24    # The single repeated root
25    r = -b / (2 * a)
26
27    # Apply initial conditions to (C1 + C2 t) e^{rt}
28    # y(0)  = C1
29    # y'(0) = C2 + r * C1
30    C1 = y0
31    C2 = yp0 - r * y0
32
33    y = (C1 + C2 * t) * np.exp(r * t)
34    return y, {"r": r, "C1": C1, "C2": C2}
35
36
37# ---- worked example: y'' + 4 y' + 4 y = 0,  y(0)=1, y'(0)=-3 ----
38t = np.linspace(0, 4, 200)
39y, info = solve_repeated_root(a=1, b=4, c=4, y0=1.0, yp0=-3.0, t=t)
40
41print(info)                # {'r': -2.0, 'C1': 1.0, 'C2': -1.0}
42print("y(0)   =", y[0])    # 1.0
43print("y(0.5) =", y[(t < 0.5).sum() - 1])
44print("zero at t =", -info["C1"] / info["C2"])   # 1.0

Numerical Check with RK4

Two analytic derivations are convincing. Three independent computer verifications are conclusive. Here we re-derive the same solution by integrating the equivalent first-order vector ODE $Y' = A Y$ with classical RK4. If the numerical and analytic answers disagree, one of them is wrong; we will find that they agree to ~1e-8 on a fine grid.

RK4 on the equivalent first-order companion system

🐍repeated_root_rk4.py

Explanation(9)

Code(47)

1Imports

NumPy for arrays and the 2×2 matrix A; matplotlib for the comparison plot at the end.

4RK4 step

Classical 4-stage Runge–Kutta. It samples the right-hand side at four carefully placed points (k1 at the start, k2 and k3 at the midpoint, k4 at the end) and combines them with weights 1, 2, 2, 1. The local truncation error is O(h⁵); the global error is O(h⁴) — orders of magnitude better than Euler for the same step size.

11Reduce 2nd-order to 1st-order

Numerical integrators consume FIRST-order systems. Introduce v = y'. Then y' = v and v' = −(b/a)v − (c/a)y. Stack into Y = [y, v]: this gives Y' = A·Y, where A is the 2×2 'companion matrix'. The repeated-root condition Δ = 0 manifests here as A having a repeated eigenvalue — a Jordan block of size 2.

17The companion matrix A

Top row [0, 1] reads 'y' = v'. Bottom row encodes the original ODE. The eigenvalues of A satisfy det(A − λI) = λ² − tr(A)λ + det(A) = λ² + (b/a)λ + (c/a) = 0 — which, after multiplying by a, is EXACTLY the characteristic equation. So 'the repeated root r' and 'the repeated eigenvalue of A' are the same number, told two different ways.

21Right-hand side f(t, Y)

Because the ODE is autonomous (no explicit t) and linear, f(t, Y) is just the matrix-vector product A @ Y. The @ symbol is NumPy's matrix-multiplication operator.

24Initial state Y(0)

The two ICs become the two components of Y₀. Without BOTH y(0) and y'(0), the problem is under-determined: many different (C₁, C₂) pairs would match a single y(0).

27Time-stepping loop

Walk the grid one step at a time, applying RK4. Each iteration only needs the previous state and the local step size h. This is how every general-purpose ODE solver works under the hood — scipy.integrate.solve_ivp, MATLAB's ode45, and torchdiffeq's odeint are all dressed-up versions of this loop.

31Return only y

We computed both y(t) and v(t) = y'(t), but the user asked for y, so we slice column 0. (For mechanics problems you usually want both — y is position, v is velocity.)

35Compare against the closed form

The analytic answer is y(t) = (1 − t)e^(−2t). RK4 with 400 steps on [0, 4] matches it to ~1e-8 — well under any plotting tolerance. This is the gold-standard sanity check: numerical and analytic methods must agree on problems where both are valid.

38 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def rk4_step(f, y, t, h):
5    """One classical Runge–Kutta-4 step for y' = f(t, y)."""
6    k1 = f(t,         y)
7    k2 = f(t + h/2,   y + h * k1 / 2)
8    k3 = f(t + h/2,   y + h * k2 / 2)
9    k4 = f(t + h,     y + h * k3)
10    return y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
11
12def simulate_2nd_order(a, b, c, y0, yp0, t):
13    """
14    Convert a*y'' + b*y' + c*y = 0 into a 1st-order system
15        Y' = A @ Y,    Y = [y, y']
16    and march RK4 forward.
17    """
18    A = np.array([[0.0,        1.0],
19                  [-c/a,   -b/a]])
20
21    def f(_t, Y):
22        return A @ Y
23
24    Y = np.array([y0, yp0], dtype=float)
25    out = np.zeros((len(t), 2))
26    out[0] = Y
27    for i in range(1, len(t)):
28        h = t[i] - t[i-1]
29        Y = rk4_step(f, Y, t[i-1], h)
30        out[i] = Y
31    return out[:, 0]   # return y only
32
33# ---- numerical vs analytic on the repeated-root example ----
34t = np.linspace(0, 4, 400)
35y_num = simulate_2nd_order(a=1, b=4, c=4, y0=1.0, yp0=-3.0, t=t)
36y_exact = (1 - t) * np.exp(-2 * t)
37
38err = np.max(np.abs(y_num - y_exact))
39print("max abs error:", err)         # ~1e-8 on a fine grid
40
41# Plot both — they should overlap.
42plt.plot(t, y_exact, "k--", label="analytic (1 − t)e^(−2t)")
43plt.plot(t, y_num,         label="RK4 on Y' = AY")
44plt.axhline(0, color="gray", lw=0.7)
45plt.legend(); plt.xlabel("t"); plt.ylabel("y(t)")
46plt.title("Critically damped:  y'' + 4y' + 4y = 0,  y(0)=1, y'(0)=-3")
47plt.show()

PyTorch View: A Jordan Block in Disguise

Now the punchline. Promote the companion matrix to a torch tensor and compute $e^{A t}$ directly. The matrix exponential is defined by the same Taylor series as the scalar one,

e^{A t} = \sum_{k=0}^\infty \frac{(A t)^k}{k!} = I + A t + \tfrac{1}{2}(A t)^2 + \cdots

For a generic matrix this series is infinite. But here $A$ has a single repeated eigenvalue $r$ . The matrix $N := A - r I$ turns out to be nilpotent of order 2: $N^2 = 0$ . Using $A = r I + N$ and $e^{r I t} = e^{r t} I$ :

e^{A t} = e^{r t}\,e^{N t} = e^{r t}\bigl(I + N t + \tfrac{1}{2}(N t)^2 + \cdots\bigr) = e^{r t}\bigl(I + N t\bigr).

The series truncates at the linear term! That single $N t$ is the structural origin of the $t$ factor in the scalar solution. Repeated eigenvalue ⟺ nilpotent $N$ ⟺ matrix exponential carries a finite polynomial in $t$ ⟺ $t e^{rt}$ appears in the solution.

Closed-form solution via the matrix exponential of a Jordan block

🐍repeated_root_matrix_exp.py

Explanation(8)

Code(39)

1Import PyTorch

We use torch.linalg.matrix_exp, which evaluates the matrix exponential to machine precision via Pade approximation + scaling-and-squaring. It is the same routine used inside neural-ODE libraries for any linear-flow component.

13Build A

Same companion matrix as the NumPy version. We use float64 so the residual against the analytic answer can be ~1e-14, not 1e-7.

15Initial vector Y(0)

Stack y0 and yp0 as a 2-vector. The whole ODE has now been reduced to one linear-algebra question: 'apply exp(A·t) to Y0.'

19The nilpotent piece N = A − r·I

Here is the deepest 'why' of the t factor. When A has a single eigenvalue r of algebraic multiplicity 2 but only one independent eigenvector, A is similar to a 2×2 Jordan block, and (A − r·I)² is the zero matrix — that is, A − r·I is NILPOTENT of order 2.

20Print N @ N

PyTorch will print a 2×2 matrix of zeros (up to round-off). This is the linear-algebra fingerprint of a repeated eigenvalue: the matrix exponential series exp(A t) = Σ (A t)^k / k! becomes a FINITE polynomial. All terms beyond the first power of N vanish.

24Vectorized matrix_exp over the time grid

For every time t_i, we form A·t_i, take its matrix exponential, and multiply by Y0. Conceptually this is just Y(t_i) = exp(A t_i) Y(0); operationally it is exact (no time integration). The list comprehension stacks the 2-vectors into an (T, 2) tensor.

28Take the y component

Column 0 of the stacked tensor is the y(t) trajectory. (Column 1 would be y'(t) — useful for phase-portrait plots.)

34Compare against (1 − t)e^(−2t)

We get the same hand-computed answer (1 − t)e^(−2t) to ~1e-14. The bigger lesson: matrix_exp + Jordan structure is a single unified picture for the WHOLE classification (distinct / repeated / complex). The 't' factor is not a magic guess — it is the t¹ term you get when you exponentiate a nilpotent operator.

31 lines without explanation

1import torch
2
3def solve_repeated_via_matrix_exp(a, b, c, y0, yp0, t):
4    """
5    Solve a*y'' + b*y' + c*y = 0 (Δ = 0) using the matrix exponential.
6
7    The vector form is  Y' = A @ Y  with  Y = [y, y'], and the exact
8    solution is  Y(t) = exp(A * t) @ Y(0).
9
10    For a repeated eigenvalue r, A - r*I is NILPOTENT of order 2, so
11        exp(A t) = e^{rt} * (I + (A - rI) * t)
12    which is exactly where the linear-in-t factor of the closed-form
13    solution comes from.
14    """
15    A = torch.tensor([[0.0, 1.0],
16                      [-c / a, -b / a]], dtype=torch.float64)
17    Y0 = torch.tensor([y0, yp0], dtype=torch.float64)
18
19    # Compute the repeated eigenvalue analytically for the sanity print
20    r = -b / (2.0 * a)
21    I = torch.eye(2, dtype=torch.float64)
22    N = A - r * I                       # Should be NILPOTENT: N @ N == 0
23    print("N @ N (should be all zeros):")
24    print(N @ N)
25
26    # Y(t_i) = exp(A * t_i) @ Y0  for every time t_i
27    Y = torch.stack([
28        torch.linalg.matrix_exp(A * ti) @ Y0
29        for ti in t
30    ])
31    return Y[:, 0], r                    # return y(t) and the repeated root
32
33t = torch.linspace(0, 4, 200, dtype=torch.float64)
34y_pt, r = solve_repeated_via_matrix_exp(1, 4, 4, 1.0, -3.0, t)
35
36y_exact = (1.0 - t) * torch.exp(-2.0 * t)
37print("repeated root r =", r)
38print("max abs error  =",
39      torch.max(torch.abs(y_pt - y_exact)).item())

Why this view scales

For higher-order linear ODEs and for linear systems with several variables, the same picture works. An eigenvalue $\lambda$ with algebraic multiplicity $m$ but only $k < m$ independent eigenvectors produces solutions of the form $(p_0 + p_1 t + \cdots + p_{m-k}\,t^{m-k})\,e^{\lambda t}$ . The same Jordan-block mechanism explains $t e^{rt}$ in 2nd order, and it explains $t^2 e^{rt}$ the first time you meet a triple root in a 3rd-order equation.

Common Pitfalls

Writing the wrong second solution

The most common mistake is to write $y = C_1 e^{rt} + C_2 e^{rt}$ in the repeated-root case. That is the same exponential twice — the two terms combine into $(C_1 + C_2) e^{rt}$ , a ONE-parameter family. The correct second solution is $t e^{rt}$ , not another $e^{rt}$ .

Forgetting the product rule for y'(0)

Differentiating $(C_1 + C_2 t) e^{rt}$ requires the product rule. A frequent error is to write $y'(0) = C_2$ , dropping the $r C_1$ term. Always carry both pieces: $y'(0) = C_2 + r C_1$ .

Mistaking Δ ≈ 0 for Δ = 0

With finite-precision arithmetic, you will rarely hit $\Delta = 0$ exactly. If your code branches on the sign of $\Delta$ , add a small tolerance and treat near-zero values as the repeated case — or better, use the matrix-exponential approach, which handles all three regimes uniformly.

Sign of r vs sign of b

$r = -b/(2a)$ . For $a > 0$ (the usual physical case), the solution decays when $b > 0$ (positive damping) and grows when $b < 0$ (negative damping — an energy-pumping system). Always sanity-check by reading off the sign of $r$ .

Summary

For a homogeneous second-order linear ODE with constant coefficients and discriminant $\Delta = 0$ :

The single repeated root is $r = -b/(2a)$ .
The general solution is $y(t) = (C_1 + C_2 t)\,e^{rt}$ . The mysterious factor of $t$ arises from three independent viewpoints: coalescence of two close roots, reduction of order, and the Jordan-block structure of the companion matrix.
The two constants are $C_1 = y(0)$ and $C_2 = y'(0) - r\,y(0)$ .
Physically this is the boundary between overdamped and underdamped behaviour — critical damping, $\zeta = 1$ .
The solution can change sign at most once (no oscillation), at $t = -C_1/C_2$ .

Key formulas at a glance

Quantity	Formula
Discriminant	Δ = b² − 4ac = 0
Repeated root	r = −b / (2a)
Fundamental set	{ e^(rt), t·e^(rt) }
General solution	y(t) = (C₁ + C₂t) e^(rt)
First constant	C₁ = y(0)
Second constant	C₂ = y'(0) − r·y(0)
Wronskian	W(t) = e^(2rt) (never zero)
Single zero crossing	t = −C₁/C₂ (if C₂ ≠ 0)
Matrix view	exp(A t) = e^(rt) (I + (A − rI) t), (A − rI)² = 0

The slogan

"A repeated root gives one exponential and one t times the same exponential — because the second solution is the derivative of e^(st) with respect to s, evaluated at the root."

Coming next: Section 22.04 turns to nonhomogeneous equations $a y'' + b y' + c y = f(t)$ and the method of undetermined coefficients. The homogeneous theory we have built — distinct, repeated, complex — provides the natural building blocks for the general nonhomogeneous solution.