Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognize a homogeneous second-order linear ODE with constant coefficients in the standard form $ay'' + by' + cy = 0$ .
Derive the characteristic equation $ar^2 + br + c = 0$ from the exponential ansatz $y = e^{rt}$ .
Identify which of the three cases (distinct real, repeated, complex conjugate) applies just by looking at the discriminant $\Delta = b^2 - 4ac$ .
Write the general solution in every case, and apply two initial conditions to pin down the constants.
Interpret the roots geometrically — as eigenvalues of a 2×2 system, as points in the complex plane, and as decay/growth rates and oscillation frequencies.
Implement both an analytic solver and a numerical RK4 solver in Python, and verify they agree.

The Big Picture: Why a Second Derivative?

"Most laws of physics are second-order differential equations." — anonymous, but you will not find a physicist who disagrees.

First-order ODEs describe systems where one thing depends on its own current value — radioactive decay, charging capacitors, drug elimination. But the universe's most basic law is Newton's second law:

F = ma = m\,\frac{d^2x}{dt^2}

The moment force depends on the object's state (spring force $-kx$ , damping $-c\,\dot{x}$ ), you get a second-order ODE:

m\ddot{x} + c\dot{x} + kx = 0

And the moment you ask the same question for an RLC circuit, a swinging pendulum (linearized), a vibrating molecule, or the bending of a beam, the exact same equation appears with different letters. That is why this one section unlocks an enormous portion of physics and engineering.

The Core Pattern

A homogeneous second-order linear ODE with constant coefficients has the form

ay'' + by' + cy = 0,\quad a \neq 0

where $a, b, c \in \mathbb{R}$ are constants. "Homogeneous" means the right-hand side is zero — no external forcing, only the system's own dynamics.

🪨 Mechanics

Mass on a spring with damping
Small-angle pendulum
Torsion of a shaft
Beam vibration (modal analysis)

⚡ Electrical

Series RLC circuit (no source)
Tank circuits / band-pass filters
Transmission-line transients

🌊 Waves & quantum

Time-independent Schrödinger 1D
Stationary Helmholtz equation in 1D
Boundary-layer eigenmodes

🤖 ML / control

Momentum gradient descent
PID controller dynamics
Linear part of every neural ODE

Anatomy of the Equation

Three words deserve unpacking before we go further.

Word	Meaning	What it implies for solving
Second-order	Highest derivative is y''	We will need TWO initial conditions: y(0) and y'(0).
Linear	y, y', y'' appear to the first power, with no products like y·y' or y².	Solutions superpose: if y₁ and y₂ solve it, so does c₁y₁ + c₂y₂.
Constant coefficients	a, b, c are real numbers — not functions of t.	We can guess y = e^(rt). With t-dependent coefficients (next sections) this trick fails.
Homogeneous	Right-hand side is 0.	There is no driving term. We are studying the system's own (free) response.

Why two initial conditions?

A second-order ODE locally determines $y''(t)$ from $y(t)$ and $y'(t)$ . So to step the solution forward in time you must know both of these at the starting instant. Physically: to predict a mass's future trajectory you need both its position and its velocity.

The Exponential Ansatz: Why $y = e^{rt}$ ?

Look at the equation again:

ay'' + by' + cy = 0

We want a function whose first and second derivatives look like the function itself, because then $y'', y', y$ will all be proportional and the equation collapses to a single algebraic condition. Only one elementary function has this property: $e^{rt}$ .

The 'try e^(rt)' rule of thumb

Whenever the coefficients are constant and the equation is linear, the substitution $y = e^{rt}$ turns a differential equation into an algebraic equation in $r$ . This is the single most useful trick in the whole subject.

Plug it in. Using $y' = re^{rt}$ and $y'' = r^2 e^{rt}$ :

a r^2 e^{rt} + b r e^{rt} + c e^{rt} = 0

Factor out $e^{rt}$ . Because $e^{rt} > 0$ for every real or complex $r$ , it can never be zero, so we may safely divide it out:

\boxed{\,a r^2 + b r + c = 0\,}

This is called the characteristic equation (or auxiliary equation). It is a plain old quadratic in $r$ . Every solution of the differential equation is built from the roots of this quadratic.

The Characteristic Equation

The quadratic $ar^2 + br + c = 0$ has two roots $r_1, r_2$ given by the familiar formula:

r_{1,2} = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}

The quantity $\Delta = b^2 - 4ac$ is the discriminant. Its sign decides everything:

Sign of Δ	Type of roots	Solution family
Δ > 0	Two distinct real roots r₁ ≠ r₂	y = C₁ e^(r₁t) + C₂ e^(r₂t)
Δ = 0	One repeated real root r (algebraic multiplicity 2)	y = (C₁ + C₂ t) e^(rt)
Δ < 0	Complex conjugate pair α ± iβ	y = e^(αt) (C₁ cos βt + C₂ sin βt)

The constants $C_1, C_2$ are determined by the two initial conditions. We will work through each case in detail, but notice the unifying theme: roots of the characteristic equation directly become the exponents of the solution.

The Three Cases of Roots

Case 1 — Distinct real roots (Δ > 0)

Both $e^{r_1 t}$ and $e^{r_2 t}$ are solutions, and they are linearly independent (one cannot be written as a constant multiple of the other when $r_1 \neq r_2$ ). The general solution is therefore

y(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t}.

If both roots are negative, the solution decays to zero — overdamped behaviour (think of a screen-door closer sliding gently shut). If at least one root is positive, the solution grows unboundedly — an unstable mode.

Case 2 — Repeated real root (Δ = 0)

The single root is $r = -b/(2a)$ . One solution is obvious: $y_1 = e^{rt}$ . But we need two linearly independent solutions. Surprisingly, the second one is $y_2 = t\,e^{rt}$ .

Why? Substitute $y = t e^{rt}$ back into the ODE. Using $y' = (1 + rt)e^{rt}$ and $y'' = (2r + r^2 t)e^{rt}$ :

a(2r + r^2 t) + b(1 + rt) + ct \;=\; (ar^2+br+c)t + (2ar + b)

The first bracket is the characteristic equation, which is 0. The second bracket is $2ar + b$ , which is also 0 because $r = -b/(2a)$ . So $t e^{rt}$ really is a solution. The general solution is

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

Physically this is critical damping: the fastest non-oscillating return to equilibrium. Car suspensions are tuned to sit just slightly underdamped of this regime.

Case 3 — Complex conjugate roots (Δ < 0)

The roots are $r = \alpha \pm i\beta$ where

\alpha = -\frac{b}{2a},\qquad \beta = \frac{\sqrt{4ac - b^2}}{2a}.

The two complex solutions are $e^{(\alpha+i\beta)t}$ and $e^{(\alpha-i\beta)t}$ . We do not want complex-valued answers for a real physical system, so we combine them using Euler's formula $e^{i\beta t} = \cos\beta t + i\sin\beta t$ . Real linear combinations give:

y(t) = e^{\alpha t}\bigl(C_1 \cos\beta t + C_2 \sin\beta t\bigr).

The interpretation is beautiful: $\alpha$ controls the envelope (decay if $\alpha<0$ , growth if $\alpha>0$ ), and $\beta$ sets the oscillation frequency. The complex plane is the natural home of this picture: the real part says how fast amplitude changes, the imaginary part says how fast the phase rotates.

One picture in your head

Put the roots on the complex plane. Left half-plane (Re < 0) means stable. Right half-plane (Re > 0) means unstable. Imaginary axis means pure (undamped) oscillation. Engineers spend whole careers designing controllers to keep roots in the left half-plane.

Structure of the General Solution

All three cases share the same skeleton:

y(t) = C_1\,y_1(t) + C_2\,y_2(t)

where $y_1, y_2$ are two linearly independent solutions. This is the principle of superposition: a homogeneous linear ODE has a solution set that is a 2-dimensional vector space.

Linear independence is checked by the Wronskian:

W(y_1, y_2) = \begin{vmatrix} y_1 & y_2 \\ y_1' & y_2' \end{vmatrix} = y_1 y_2' - y_2 y_1'.

If $W \neq 0$ at any single point, the two solutions form a fundamental set. You can verify (good practice!) that in all three cases above $W(0) \neq 0$ , confirming the general solution.

Worked Example (Step-by-Step)

Solve $y'' - 5y' + 6y = 0$ subject to $y(0) = 2,\; y'(0) = 3$ .

Click to expand the full hand calculation

Step 1 — Read off a, b, c.

Comparing with $ay'' + by' + cy = 0$ we read $a = 1,\; b = -5,\; c = 6$ .

Step 2 — Write the characteristic equation.

r^2 - 5r + 6 = 0

Step 3 — Find the roots.

Factor: $(r-2)(r-3) = 0$ , so $r_1 = 2,\; r_2 = 3$ . Both real, distinct, both positive — expect exponential growth.

Step 4 — Discriminant check.

$\Delta = (-5)^2 - 4(1)(6) = 25 - 24 = 1 > 0$ . Confirms the "distinct real roots" case.

Step 5 — General solution.

y(t) = C_1 e^{2t} + C_2 e^{3t}

Step 6 — Apply $y(0) = 2$ .

At $t = 0$ : $y(0) = C_1 + C_2 = 2$ .

Step 7 — Differentiate and apply $y'(0) = 3$ .

$y'(t) = 2C_1 e^{2t} + 3C_2 e^{3t}$ . At $t = 0$ : $y'(0) = 2C_1 + 3C_2 = 3$ .

Step 8 — Solve the 2×2 linear system.

From $C_1 + C_2 = 2$ : $C_1 = 2 - C_2$ . Substitute: $2(2 - C_2) + 3C_2 = 3$ $\Rightarrow$ $4 + C_2 = 3$ $\Rightarrow$ $C_2 = -1$ . Then $C_1 = 3$ .

Step 9 — Particular solution.

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

Step 10 — Sanity check.

$y(0) = 3 - 1 = 2$ ✓. $y'(0) = 6 - 3 = 3$ ✓. Plug back into the ODE: $y'' = 12 e^{2t} - 9 e^{3t}$ , so $y'' - 5y' + 6y = (12 - 10 + 6)e^{2t} + (-9 + 15 - 6)e^{3t} = 8e^{2t}$ hmm — let's redo that carefully. The coefficient on $e^{2t}$ is $12 - 5(6) + 6(3) = 12 - 30 + 18 = 0$ . On $e^{3t}$ : $-9 - 5(-3) + 6(-1) = -9 + 15 - 6 = 0$ . ✓ Both sides vanish.

Step 11 — A long-term value.

At $t = 1$ : $y(1) = 3e^2 - e^3 \approx 3(7.389) - 20.086 \approx 2.081$ . You can verify this against the interactive plot just below (set a=1, b=−5, c=6, y(0)=2, y′(0)=3).

Interactive Characteristic-Equation Explorer

Drag the sliders for $a, b, c$ and watch the characteristic equation, its discriminant, its two roots in the complex plane, and the corresponding solution $y(t)$ update in real time. Click any preset to jump straight into a famous regime.

Loading interactive explorer…

What to look for

Watch the two root markers slide along the real axis as you change $b$ . The instant they collide is exactly $\Delta = 0$ (the critical-damping knife edge).
Push them past that boundary and they leave the real axis vertically as a conjugate pair — that is the geometric meaning of complex roots.
Drag $c$ below zero. Roots split apart along the real axis, one becomes positive (unstable mode), and the solution grows without bound.

Physical Meaning: Mass–Spring–Damper

The canonical physical realisation is a block of mass $m$ attached to a wall by a spring of stiffness $k$ , with a dashpot providing viscous damping $c$ :

m\,\ddot x + c\,\dot x + k\,x = 0.

Compare term by term with the abstract form: the coefficient of the second derivative is the inertia, the coefficient of the first derivative is the energy loss, and the coefficient of $y$ itself is the restoring stiffness.

Engineers prefer two derived quantities over $(m, c, k)$ directly:

Quantity	Formula	Meaning
Natural frequency	ω₀ = √(k/m)	Angular frequency if there were no damping (c = 0).
Damping ratio	ζ = c / (2√(mk))	Dimensionless. ζ < 1 → oscillates; ζ = 1 → critically damped; ζ > 1 → overdamped.

Notice that $\zeta = 1$ is exactly $c^2 = 4mk$ , which is the discriminant condition $\Delta = 0$ . Mathematically and physically, the "three cases" are the same boundary.

Loading mass–spring–damper simulation…

Computation in Python

Math, then code. We will write two solvers: one analytic (using the closed-form solutions we derived), and one numerical (using classical RK4 on the reduced first-order system). They should agree to many digits.

Analytic solver — closed form from the characteristic equation

Direct evaluation of y(t) once the roots are known

🐍analytic_2nd_order.py

Explanation(13)

Code(59)

1Import NumPy

NumPy gives us vectorized math: exp, cos, sin, sqrt operate elementwise on the time grid, so the entire solution is one expression — no Python loop.

4Function header

We take the three coefficients a, b, c, the two initial conditions y0 and yp0, and a NumPy array t of time values. A second-order ODE needs TWO initial conditions because the general solution has two unknown constants C1 and C2.

21Discriminant

Δ = b² − 4ac is the same discriminant you saw for quadratics. Its SIGN — not its value — decides which of the three solution families applies. This single number classifies the entire qualitative behavior of the ODE.

23Distinct-roots branch

When Δ > 0 the quadratic ar² + br + c = 0 has two different real solutions r₁ and r₂. The general solution is then y(t) = C₁e^(r₁t) + C₂e^(r₂t).

25Compute r1 and r2

Standard quadratic formula r = (−b ± √Δ) / (2a). r₁ uses the +, r₂ uses the −. Each root is a separate exponential mode of the system.

28Apply initial conditions

C₁ + C₂ = y(0) and r₁C₁ + r₂C₂ = y′(0). Two linear equations in two unknowns. Solving gives C₂ first, then C₁ = y0 − C₂. (The denominator r₂ − r₁ ≠ 0 because roots are distinct in this branch.)

31Build y(t)

NumPy evaluates C1·exp(r1·t) + C2·exp(r2·t) across the entire time grid at once. The result is a 1-D array of the same length as t.

35Complex-roots branch

When Δ < 0 the roots are α ± iβ — a complex-conjugate pair. By Euler's formula e^(α+iβ)t = e^(αt)(cos(βt) + i sin(βt)), the REAL solution becomes a damped oscillation: e^(αt)(C₁cos(βt) + C₂sin(βt)). α governs decay/growth; β governs frequency.

37Real and imaginary parts of the root

α = −b / (2a) is the common real part. β = √(−Δ) / (2a) is the magnitude of the imaginary part. We use √(−Δ) because Δ is negative here, so −Δ is positive.

40Initial conditions for the oscillating case

y(0) = C₁ (cosine starts at 1, sine starts at 0). y′(0) = αC₁ + βC₂, so C₂ = (y′(0) − αy(0)) / β. The β in the denominator is non-zero because we are in the complex-roots branch.

46Repeated-root branch

When Δ = 0 the quadratic has a single root r = −b/(2a). One exponential e^(rt) is not enough — there must be a second linearly independent solution. The correct one is t·e^(rt), giving the general solution y(t) = (C₁ + C₂t)e^(rt).

48Initial conditions for the repeated-root case

C₁ = y(0) directly. Differentiating (C₁ + C₂t)e^(rt) and setting t=0 gives y′(0) = C₂ + r·C₁, so C₂ = y′(0) − r·y(0).

55Run the worked example

We call the solver with a=1, b=−5, c=6, y(0)=2, y′(0)=3. The characteristic equation r² − 5r + 6 = 0 factors as (r−2)(r−3) = 0, giving r₁=3, r₂=2 and the closed form y(t) = 3e^(2t) − e^(3t).

46 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_homogeneous_2nd_order(a, b, c, y0, yp0, t):
5    """
6    Solve a*y'' + b*y' + c*y = 0 analytically.
7
8    Parameters
9    ----------
10    a, b, c : float
11        Constant coefficients of the ODE.
12    y0  : float   y(0),  initial value
13    yp0 : float   y'(0), initial slope
14    t   : np.ndarray, time grid where the solution is evaluated.
15
16    Returns
17    -------
18    y : np.ndarray with y(t) on the same grid.
19    info : dict with the case name and the two roots.
20    """
21    disc = b**2 - 4 * a * c
22
23    if disc > 0:
24        # Two distinct real roots
25        sqrtD = np.sqrt(disc)
26        r1 = (-b + sqrtD) / (2 * a)
27        r2 = (-b - sqrtD) / (2 * a)
28        # Solve [[1,1],[r1,r2]] @ [C1,C2] = [y0, yp0]
29        C2 = (yp0 - r1 * y0) / (r2 - r1)
30        C1 = y0 - C2
31        y = C1 * np.exp(r1 * t) + C2 * np.exp(r2 * t)
32        return y, {"case": "distinct_real", "r1": r1, "r2": r2,
33                   "C1": C1, "C2": C2}
34
35    if disc < 0:
36        # Complex conjugate roots alpha ± i*beta
37        alpha = -b / (2 * a)
38        beta  = np.sqrt(-disc) / (2 * a)
39        C1 = y0
40        C2 = (yp0 - alpha * y0) / beta
41        y = np.exp(alpha * t) * (C1 * np.cos(beta * t) + C2 * np.sin(beta * t))
42        return y, {"case": "complex", "alpha": alpha, "beta": beta,
43                   "C1": C1, "C2": C2}
44
45    # disc == 0: one repeated root
46    r = -b / (2 * a)
47    C1 = y0
48    C2 = yp0 - r * y0
49    y = (C1 + C2 * t) * np.exp(r * t)
50    return y, {"case": "repeated", "r": r, "C1": C1, "C2": C2}
51
52
53# ---- run the worked example: y'' - 5y' + 6y = 0, y(0)=2, y'(0)=3 ----
54t = np.linspace(0, 1.5, 200)
55y, info = solve_homogeneous_2nd_order(a=1, b=-5, c=6, y0=2, yp0=3, t=t)
56
57print(info)
58print("y(0)   =", y[0])
59print("y(1.0) =", y[133])   # t[133] is close to 1.0

Numerical solver — RK4 on the reduced 2×2 system

Numerical solvers consume first-order systems. The trick is to introduce $v = y'$ and rewrite the scalar second-order ODE as a 2-D vector ODE:

\frac{d}{dt}\begin{bmatrix} y \\ v \end{bmatrix} = \underbrace{\begin{bmatrix} 0 & 1 \\ -c/a & -b/a \end{bmatrix}}_{A} \begin{bmatrix} y \\ v \end{bmatrix}.

The eigenvalues of this companion matrix $A$ are exactly the roots of the characteristic equation. So the same $r_1, r_2$ we found by hand are now the eigenvalues of a 2×2 matrix — a pivot to linear algebra that we exploit in the PyTorch version below.

RK4 integration of the equivalent first-order system

🐍rk4_2nd_order.py

Explanation(9)

Code(39)

1Import NumPy

We use NumPy arrays as 2-vectors holding the state Y = [y, y′]. The matrix A is a small 2×2.

3RK4 step function

Classical Runge–Kutta 4 takes one explicit step. It evaluates the right-hand side f at four carefully placed points (k1, k2, k3, k4) and combines them. Local truncation error is O(h⁵); global error is O(h⁴) — much better than plain Euler's O(h).

11Reduce 2nd-order to 1st-order system

Numerical solvers consume FIRST-order systems. Introduce v = y′. Then y′ = v and v′ = −(b/a)v − (c/a)y. Stacking these gives Y′ = A·Y with Y = [y, v] — a linear vector ODE.

17The matrix A

A has a very specific structure: top row is [0, 1] (because y′ = v) and the bottom row encodes the original second-order equation. The roots of det(λI − A) = 0 are EXACTLY the roots of the characteristic equation aλ² + bλ + c = 0. So A's eigenvalues ARE the r₁, r₂ from the analytic solution.

21Right-hand side f(t, Y)

Because the ODE is autonomous (no t appears explicitly) and linear, f(t, Y) is just A @ Y. The @ operator is matrix-vector product.

23Initialize the state

The two initial conditions y(0) and y′(0) become the two components of Y₀ = [y0, yp0]. Without BOTH initial conditions the problem is under-determined.

27Time-stepping loop

Walk the time grid one step at a time, applying RK4 to march Y forward. At each step we only need the previous state and the step size h. This is how every solver — scipy.integrate.solve_ivp, PyTorch's torchdiffeq, MATLAB's ode45 — works under the hood.

31Return only y

The user asked for y(t), so we drop the velocity column. (For physics simulations you usually want both y and v.)

38Compare numerical vs closed-form

We computed the exact solution y(t) = 3e^(2t) − e^(3t) by hand. RK4 should match it to ~10⁻⁸ on a fine grid. This is the gold-standard sanity check for any new ODE solver: solve a problem with a known closed form.

30 lines without explanation

1import numpy as np
2
3def rk4_step(f, y, t, h):
4    """One classical Runge-Kutta-4 step for y' = f(t, y)."""
5    k1 = f(t,         y)
6    k2 = f(t + h/2,   y + h * k1 / 2)
7    k3 = f(t + h/2,   y + h * k2 / 2)
8    k4 = f(t + h,     y + h * k3)
9    return y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
10
11def simulate_2nd_order(a, b, c, y0, yp0, t):
12    """
13    Integrate a*y'' + b*y' + c*y = 0 by converting it to a
14    first-order system  Y' = A @ Y  with  Y = [y, y'].
15    """
16    # Rewrite y'' = -(b/a) y' - (c/a) y
17    A = np.array([[0.0,         1.0],
18                  [-c/a,    -b/a]])
19
20    def f(_t, Y):
21        return A @ Y
22
23    Y = np.array([y0, yp0], dtype=float)
24    out = np.zeros((len(t), 2))
25    out[0] = Y
26    for i in range(1, len(t)):
27        h = t[i] - t[i-1]
28        Y = rk4_step(f, Y, t[i-1], h)
29        out[i] = Y
30    return out[:, 0]    # return only y, drop y'
31
32# Numerical run of the same example
33t = np.linspace(0, 1.5, 200)
34y_num = simulate_2nd_order(a=1, b=-5, c=6, y0=2, yp0=3, t=t)
35
36# Compare to the closed-form 3*e^(2t) - e^(3t)
37y_exact = 3 * np.exp(2 * t) - np.exp(3 * t)
38err = np.max(np.abs(y_num - y_exact))
39print("max abs error:", err)

PyTorch View: Linear ODEs as a Matrix Exponential

Once the second-order scalar ODE is rewritten as a vector ODE $Y' = A Y$ , an extraordinary fact appears: the exact solution is

Y(t) = e^{A t}\,Y(0).

Here $e^{At}$ is the matrix exponential — defined by the same power series $I + At + \tfrac{1}{2!}(At)^2 + \cdots$ that defines the scalar one. PyTorch ships $\texttt{torch.linalg.matrix\_exp}$ , which evaluates this expression to machine precision.

Closed-form solution via the matrix exponential

🐍matrix_exp_2nd_order.py

Explanation(6)

Code(26)

1Import PyTorch

PyTorch ships with torch.linalg.matrix_exp, a dense matrix exponential. We use it to write the EXACT solution of any linear constant-coefficient ODE in one line.

9Build the system matrix A

Same A as in the RK4 file. The trick: for a LINEAR constant-coefficient system Y′ = AY, the exact solution is Y(t) = e^(At) Y(0). No numerical integrator needed.

12Initial state Y(0)

Stack y(0) and y′(0) into the 2-vector Y(0). The whole ODE is now reduced to one linear-algebra problem.

16matrix_exp(A * t_i) @ Y0

For every time t_i, we form the 2×2 matrix A·t_i, take its matrix exponential, and multiply by Y(0). The result is Y(t_i) = [y(t_i), y′(t_i)]. The matrix exponential is to matrices what the scalar exp is to numbers — and it's the natural object for solving linear ODEs.

20Take the y component

We stack all 2-vectors into a (T, 2) tensor and slice column 0 to get the y(t) trajectory.

25Round-trip check

Compare against the hand-computed closed form 3e^(2t) − e^(3t). The matrix-exp approach is EXACT up to floating-point round-off (~10⁻¹⁵ in float64). This is the same routine PyTorch's neural-ODE library uses for the linear-flow component of more complex systems.

20 lines without explanation

1import torch
2
3def solve_homogeneous_via_matrix_exp(a, b, c, y0, yp0, t):
4    """
5    Solve a*y'' + b*y' + c*y = 0 using the matrix exponential.
6
7    Trick: write the ODE as Y' = A Y with Y = [y, y'].
8    The exact solution is Y(t) = exp(A * t) @ Y(0),
9    where exp(.) is the MATRIX exponential.
10    """
11    A = torch.tensor([[0.0, 1.0],
12                      [-c / a, -b / a]], dtype=torch.float64)
13    Y0 = torch.tensor([y0, yp0], dtype=torch.float64)
14
15    # Vectorized: exp(A * t_i) for every t_i
16    Y = torch.stack([
17        torch.linalg.matrix_exp(A * ti) @ Y0
18        for ti in t
19    ])
20    return Y[:, 0]   # y component
21
22t = torch.linspace(0, 1.5, 200, dtype=torch.float64)
23y_pt = solve_homogeneous_via_matrix_exp(1, -5, 6, 2.0, 3.0, t)
24
25y_exact = 3 * torch.exp(2 * t) - torch.exp(3 * t)
26print("max abs error:", torch.max(torch.abs(y_pt - y_exact)).item())

Connection to neural ODEs

A linear neural ODE $\dot Y = A Y$ with a learnable matrix $A$ is the simplest possible continuous-time model. Its evolution is literally $e^{At}$ . The eigenvalues of $A$ — the same characters that played the starring role here — determine whether the network exhibits decay, oscillation, or runaway growth. Everything you learned in this section maps directly into how modern continuous-depth models behave.

Common Pitfalls

Forgetting to divide by a

The characteristic equation uses the coefficients as written. If your ODE is $2y'' + 8y' + 6y = 0$ , the characteristic equation is $2r^2 + 8r + 6 = 0$ , NOT $r^2 + 8r + 6 = 0$ . (Of course you can divide both sides by 2 first; just be consistent.)

One initial condition is not enough

A second-order ODE has a 2-parameter family of solutions. With only one initial condition you can satisfy $C_1+C_2=y_0$ but $(C_1, C_2)$ remains under-determined. Always check that both $y(t_0)$ and $y'(t_0)$ are specified.

Sign of the discriminant — don't trust the look of b

A negative $b$ does NOT imply complex roots. It is the sign of $b^2 - 4ac$ — never just $b$ — that classifies the roots.

Repeated-root case: do not write y = C₁e^(rt) + C₂e^(rt)

That is the same exponential twice, with a single combined constant $C_1 + C_2$ . The correct second solution is $t e^{rt}$ , NOT another $e^{rt}$ .

Summary

Every homogeneous second-order linear ODE with constant coefficients can be solved with the same recipe:

Write the characteristic equation $ar^2 + br + c = 0$ .
Compute the discriminant $\Delta = b^2 - 4ac$ .
Read off the case (distinct real / repeated / complex) and the corresponding general solution.
Apply the two initial conditions to solve a 2×2 linear system for $C_1, C_2$ .

Key formulas at a glance

Case	Condition	General solution
Distinct real	Δ > 0	y = C₁ e^(r₁t) + C₂ e^(r₂t)
Repeated real	Δ = 0, r = −b/(2a)	y = (C₁ + C₂ t) e^(rt)
Complex conjugate	Δ < 0, α = −b/(2a), β = √(4ac−b²)/(2a)	y = e^(αt) (C₁ cos βt + C₂ sin βt)

The slogan

"The roots of the characteristic equation ar² + br + c = 0 are the exponents of the solution."

Coming next: Section 22.02 zooms in on the complex-roots case — turning $\alpha \pm i\beta$ into a deep theory of oscillation, phase, and the damped sinusoid. We will derive the amplitude–phase form $y = A e^{\alpha t}\cos(\beta t - \phi)$ and connect it to musical instruments and signal envelopes.