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 homogeneous linear ODE has complex conjugate roots just by checking the sign of $\Delta = b^2 - 4ac$ .
Derive the real-valued general solution $y(t) = e^{\alpha t}\bigl(C_1\cos(\beta t) + C_2\sin(\beta t)\bigr)$ from the complex form using Euler's formula.
Interpret the real part $\alpha$ as decay/growth rate and the imaginary part $\beta$ as angular frequency.
Rewrite the solution in amplitude-phase form $y(t) = A\,e^{\alpha t}\cos(\beta t - \varphi)$ and read off period, half-life, and phase.
Translate between coefficient form $(a,b,c)$ and physical form $(\omega_n,\zeta)$ , then classify damping regimes.
Implement a complex-root solver in plain Python and cross-check it against PyTorch's $\mathrm{matrix\_exp}$ .

Where We Pick Up: A New Kind of Root

In section 22.1 we tried $y = e^{rt}$ in $ay'' + by' + cy = 0$ and ended up with the characteristic equation:

a r^2 + b r + c = 0

The quadratic formula gave us three branches depending on $\Delta = b^2 - 4ac$ . We handled the two easy cases — distinct real roots and the repeated real root. This section is the third and richest one:

\Delta < 0 \quad\Longrightarrow\quad r = \alpha \pm i\beta,\ \beta > 0.

The promise: negative discriminants do not break the method — they turn the boring exponential growth/decay into something far more interesting: oscillation. Every spring, every pendulum, every RLC circuit, every vibrating molecule lives here.

Two questions are about to dominate the whole section:

We started in $\mathbb{R}$ . How did $i = \sqrt{-1}$ sneak in, and what do we do with it?
Why does an imaginary root manifest as a real oscillation? What is the bridge?

The bridge is Euler's formula. Once it is in place, the rest of the section is just careful bookkeeping.

The Imaginary Wall

Take the cleanest possible example: $y'' + y = 0$ . The characteristic equation is $r^2 + 1 = 0$ , giving $r = \pm i$ .

Naively plug into the "general solution" from section 22.1:

y(t) = C_1 e^{i t} + C_2 e^{-i t}.

This looks complex-valued. But we are modelling a real physical thing — a frictionless mass on a spring, or an LC circuit. Voltages and positions are real numbers. We need a way to express $e^{it}$ in real terms.

The discomfort is the point

The complex form is not wrong. It is simply a coordinate system that is awkward for a real-valued physical observable. The job of this section is to translate the complex coordinates into a familiar real basis — sines and cosines — without losing any information.

Euler's Bridge — From $e^{i\beta t}$ to Sines and Cosines

The single most important identity in this whole chapter is Euler's formula:

e^{i\theta} = \cos\theta + i\sin\theta.

Why does this hold?

Expand the Taylor series. The series for $e^{i\theta}$ can be regrouped into the series for $\cos\theta$ and $i\sin\theta$ :

e^{i\theta} = \sum_{n=0}^{\infty} \frac{(i\theta)^n}{n!} = \underbrace{1 - \frac{\theta^2}{2!} + \frac{\theta^4}{4!} - \cdots}_{\cos\theta} + i\underbrace{\Bigl(\theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots\Bigr)}_{\sin\theta}

The trick is that $i^2 = -1$ , so even powers of $i$ alternate sign — exactly matching the alternating signs in the cosine and sine series.

Apply it to $\theta = \beta t$ :

e^{i\beta t} = \cos(\beta t) + i\sin(\beta t),\qquad e^{-i\beta t} = \cos(\beta t) - i\sin(\beta t).

Now stitch the decaying envelope back in. For a root $r = \alpha + i\beta$ :

e^{rt} = e^{(\alpha + i\beta)t} = e^{\alpha t}\bigl(\cos(\beta t) + i\sin(\beta t)\bigr).

The exponential $e^{rt}$ has two parts: a real exponential $e^{\alpha t}$ that decays ( $\alpha < 0$ ) or grows ( $\alpha > 0$ ), and a complex rotation $e^{i\beta t}$ that spins around the unit circle at angular speed $\beta$ . Multiplied together they trace out a spiral in the complex plane.

The Real-Valued General Solution

Our complex pair is $r_1 = \alpha + i\beta$ , $r_2 = \alpha - i\beta$ . By linearity,

y(t) = K_1 e^{(\alpha + i\beta)t} + K_2 e^{(\alpha - i\beta)t}

solves the ODE for any complex $K_1, K_2$ . To make $y(t)$ real for all $t$ , we choose $K_2 = \overline{K_1}$ (complex conjugate). With $K_1 = \tfrac{1}{2}(C_1 - iC_2)$ and $K_2 = \tfrac{1}{2}(C_1 + iC_2)$ , the imaginary parts cancel and you are left with this clean real form:

General solution — complex-roots case

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

where $\alpha = -\tfrac{b}{2a}$ , $\beta = \tfrac{\sqrt{-\Delta}}{2a}$ , and $C_1, C_2 \in \mathbb{R}$ are determined by the initial conditions.

Reading the formula

The $e^{\alpha t}$ in front is the envelope. Whatever happens inside the parentheses, its amplitude is bounded by $A\,e^{\alpha t}$ . The inside is a weighted sum of a cosine and a sine, both with the SAME angular frequency $\beta$ .

Why two real constants, not four?

A second-order ODE has a two-dimensional solution space. $\cos(\beta t)$ and $\sin(\beta t)$ are linearly independent and span that space (after multiplying by the envelope). Two constants $C_1, C_2$ are exactly enough to match the two initial conditions $y(0)$ and $y'(0)$ .

Applying the initial conditions

At $t = 0$ :

$y(0) = e^{0}(C_1\cdot 1 + C_2\cdot 0) = C_1$ , so $C_1 = y(0)$ .
$y'(t) = e^{\alpha t}\bigl[(\alpha C_1 + \beta C_2)\cos(\beta t) + (\alpha C_2 - \beta C_1)\sin(\beta t)\bigr]$ , so $y'(0) = \alpha C_1 + \beta C_2$ .

Solving for $C_2$ :

C_1 = y(0),\qquad C_2 = \frac{y'(0) - \alpha\,y(0)}{\beta}.

The denominator $\beta$ is non-zero exactly because we are in the complex-roots branch.

Anatomy of $\alpha$ and $\beta$

These two numbers are everything. They live in different directions of the complex plane and they tell different physical stories.

Quantity	Formula	Physical meaning	Geometric meaning
α (real part)	−b / (2a)	Decay rate of the envelope. α < 0 ⇒ amplitude shrinks; α > 0 ⇒ amplitude grows.	Horizontal position of the root in the complex plane. Left half-plane = stable.
β (imag part)	√(−Δ) / (2a)	Angular frequency in rad/s. Period T = 2π/β, frequency f = β/(2π).	Distance of the root from the real axis. Far from the axis = fast wiggle.
α² + β²	c / a	Distance² of the root from the origin (in the complex plane).	For a spring this equals ωₙ², the natural frequency squared.
2α	−b / a	Sum of the two roots (Vieta's formula).	Twice the real part — used as the trace of the system matrix.

A two-knob mental model

Imagine you are designing an instrument. The knob $\alpha$ controls how quickly the note fades; the knob $\beta$ controls the pitch. Set $\alpha = 0$ and the note rings forever. Make $\alpha$ more negative and the note dies sooner. Change $\beta$ and the same envelope holds a different musical pitch.

Amplitude–Phase Form: One Tilted Cosine

The form $C_1\cos(\beta t) + C_2\sin(\beta t)$ is algebraically convenient but visually awkward — it's a SUM of two wiggles. There is a beautiful alternative.

Use the identity $C_1\cos\theta + C_2\sin\theta = A\cos(\theta - \varphi)$ with

A = \sqrt{C_1^2 + C_2^2},\qquad \varphi = \mathrm{atan2}(C_2,\, C_1).

Then the whole solution collapses to a single damped cosine:

\boxed{\,y(t) = A\,e^{\alpha t}\cos(\beta t - \varphi)\,}

Reading this off a measured oscilloscope trace is now trivial:

Amplitude $A$ = the initial peak (before any decay).
Envelope $A\,e^{\alpha t}$ = the dashed curve touching every peak.
Frequency $\beta$ = how fast it wiggles (rad/s).
Phase $\varphi$ = how far the first peak is shifted from $t = 0$ .

Why atan2, not atan?

$\mathrm{atan2}(y, x)$ uses the signs of both arguments to return an angle in $(-\pi, \pi]$ , distinguishing all four quadrants. Plain $\arctan(C_2/C_1)$ would collapse the second and fourth quadrants together (since dividing two negatives kills the sign info), giving you the wrong phase half the time.

Worked Example (Step-by-Step)

Take this concrete IVP:

y'' + 2y' + 5y = 0,\qquad y(0) = 1,\quad y'(0) = 0.

We will derive every number on paper, then verify with the explorer and the Python solver below.

▸ Show full hand-derivation (recommended on first read)

Step 1 — Characteristic equation

Substitute $y = e^{rt}$ :

r^2 + 2r + 5 = 0.

Step 2 — Discriminant

$\Delta = 2^2 - 4(1)(5) = 4 - 20 = -16 < 0.$ So we are in the complex-roots branch.

Step 3 — Roots

$r = \dfrac{-2 \pm \sqrt{-16}}{2} = \dfrac{-2 \pm 4i}{2} = -1 \pm 2i.$

So $\alpha = -1$ and $\beta = 2$ .

Step 4 — General solution

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

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

Plug $t = 0$ : $y(0) = 1\cdot(C_1\cdot 1 + C_2\cdot 0) = C_1$ . So $C_1 = 1$ .

Step 6 — Apply $y'(0) = 0$

From $C_2 = (y'(0) - \alpha\,y(0))/\beta$ :

C_2 = \dfrac{0 - (-1)(1)}{2} = \dfrac{1}{2}.

Step 7 — Final closed form

y(t) = e^{-t}\bigl(\cos 2t + \tfrac{1}{2}\sin 2t\bigr).

Step 8 — Amplitude-phase form

$A = \sqrt{1^2 + (1/2)^2} = \sqrt{5}/2 \approx 1.1180.$
$\varphi = \mathrm{atan2}(1/2,\, 1) \approx 0.4636\,\text{rad} \approx 26.57^\circ.$

y(t) = \tfrac{\sqrt{5}}{2}\,e^{-t}\cos(2t - 0.4636).

Step 9 — Spot-check at four times

t	y(t) by hand	envelope ±A·e^(αt)
0	1.000000	±1.118034
π/4 ≈ 0.7854	e^(−π/4)·(0 + 0.5) ≈ 0.2280	±0.5098
π/2 ≈ 1.5708	e^(−π/2)·(−1 + 0) ≈ −0.2079	±0.2324
π ≈ 3.1416	e^(−π)·(1 + 0) ≈ 0.04321	±0.04832

Every value of $y(t)$ lies inside its envelope, as it must.

Step 10 — Read off the physics

Period: $T = 2\pi/\beta = \pi \approx 3.1416\,\text{s}$ .
Half-life of the envelope: $t_{1/2} = \ln 2 / |\alpha| = \ln 2 \approx 0.693\,\text{s}$ .
How many oscillations before the envelope falls to 1% of its initial value? $t = \ln 100 \approx 4.605$ ; in that time the system completes $4.605/\pi \approx 1.47$ cycles. The response dies out fast.

Interactive: Roots in the Complex Plane

Drag the upper red root anywhere in the complex plane and watch what happens to the time-domain plot and the phase-plane spiral. The envelope $\pm A\,e^{\alpha t}$ is drawn in dashed amber so you can see exactly which curve the oscillation lives inside.

Loading complex-roots explorer…

Three experiments worth doing

Pull the root all the way to the imaginary axis ( $\alpha = 0$ ). The envelope flattens to constant ±A, the phase-plane orbit closes into an ellipse. This is the undamped regime — perpetual motion.
Cross into the right half-plane ( $\alpha > 0$ ). The envelope flips to a growing exponential, the spiral flies outward. This is an unstable system — think inverted pendulum, or a poorly tuned PID loop that oscillates more violently with every cycle.
Increase $\beta$ with $\alpha$ fixed. The envelope is unchanged but the wiggles get faster. The spiral makes more turns per inward step — high-frequency, low-damping.

Damping in Disguise: $\zeta$ and $\omega_n$

Engineers rarely speak in terms of $(a, b, c)$ . They use two physically meaningful parameters: the natural frequency $\omega_n$ and the damping ratio $\zeta$ . They rewrite the ODE in the canonical form

y'' + 2\zeta\,\omega_n\,y' + \omega_n^2\,y = 0,

which matches $ay'' + by' + cy = 0$ with $\omega_n = \sqrt{c/a}$ and $\zeta = b / (2\sqrt{ac})$ .

Computing the discriminant in this notation gives $\Delta = 4\omega_n^2(\zeta^2 - 1)$ . So the three cases reorganise themselves into a single dimensionless story:

ζ value	Δ sign	Regime	α (decay rate)	β (oscillation freq)
ζ = 0	< 0	Undamped	0	ωₙ
0 < ζ < 1	< 0	Underdamped (complex roots)	−ζ·ωₙ	ωₙ·√(1 − ζ²) ≡ ωd
ζ = 1	= 0	Critically damped	−ωₙ	0 (no oscillation)
ζ > 1	> 0	Overdamped	two real roots	0 (no oscillation)

The damped frequency ωd

When you measure a real oscillation you see $\omega_d = \omega_n\sqrt{1 - \zeta^2}$ , NOT the natural frequency $\omega_n$ . Damping slows the wiggle. For $\zeta = 0.1$ the difference is only 0.5%; for $\zeta = 0.5$ it's about 13%; at $\zeta = 1$ the wiggle disappears entirely.

Engineer's rule of thumb

$\zeta \approx 0.7$ is the industry-favourite damping ratio for control systems. It gives a fast response (small $\zeta$ means snappy) with only a single tiny overshoot (large $\zeta$ means well-behaved). The trade-off between speed and overshoot is exactly what $\zeta$ parameterises.

Interactive: Four Regimes from One Spring

Same mass-spring system, four different dampers. The four ghost traces stay visible so you can compare them against the active curve. Move $\zeta$ continuously through 0, 1, and above 1 and watch the regime change in real time.

Loading damping comparator…

Python: Closed-Form Solver from Scratch

We will build a single function that takes $(a, b, c, y_0, y'_0, t)$ and returns $y(t)$ along with the physical diagnostics — $\alpha, \beta, A, \varphi$ , period, half-life. Plain NumPy, no special-case libraries.

closed-form complex-roots solver

🐍python

Explanation(13)

Code(51)

1Import NumPy

NumPy lets us evaluate exp, cos, sin, sqrt on a whole time grid in one expression. Every formula below is written so it works on a single scalar t OR a NumPy array of times.

4Function signature

We accept the three coefficients a, b, c, the two initial conditions y0 = y(0) and yp0 = y'(0), and a NumPy array t of time points where we want the solution. Two initial conditions are required because a second-order ODE has two free constants C1, C2 in its general solution.

11Compute the discriminant

Δ = b² − 4ac is exactly the same discriminant you saw in section 22.1. Its sign decides which solution family applies. Here we assert Δ < 0 so we only use the complex-roots branch.

15Real part α of the root

α = −b / (2a) is the common real part of the conjugate pair r = α ± iβ. Physically α is the EXPONENTIAL DECAY RATE of the oscillation's envelope. Negative α ⇒ amplitude shrinks; positive α ⇒ amplitude grows.

16Imaginary part β of the root

β = √(−Δ) / (2a) is the magnitude of the imaginary part. β is the ANGULAR FREQUENCY of the oscillation in radians per second. We take √(−Δ) because Δ is negative in this branch, so −Δ is positive.

23Apply initial conditions — pin C1

Plug t=0 into y(t) = e^(αt)(C₁cos(βt) + C₂sin(βt)). At t=0, e^0=1, cos(0)=1, sin(0)=0, so y(0) = C₁ directly. No algebra needed.

24Apply initial conditions — pin C2

Differentiate y and plug in t=0. You get y'(0) = αC₁ + βC₂. Solving for C₂ gives (y'(0) − αy(0))/β. The β in the denominator is non-zero because we are in the complex-roots branch (β > 0).

27Build the real solution

y(t) = e^(αt)(C₁cos(βt) + C₂sin(βt)). NumPy evaluates this whole expression elementwise on the time grid. The result is a 1-D array of length len(t). The e^(αt) factor is the envelope; the cos+sin combination is the oscillation inside it.

30Compute the amplitude A

A = √(C₁² + C₂²) is the radius of the cosine–sine combination. It is the peak height of the dimensionless oscillation (before the envelope is multiplied in). This single number summarizes the initial 'energy' of the response.

31Compute the phase φ

φ = atan2(C₂, C₁) shifts the cosine. With A and φ the solution can be rewritten as y(t) = A·e^(αt)·cos(βt − φ), which is a SINGLE damped cosine. atan2 (not atan) is used to land in the correct quadrant regardless of signs.

35Return diagnostics

We return α, β, the constants, the amplitude, the phase, the period 2π/β, and the half-life ln(2)/(−α) (for a stable decay). These quantities are what an engineer actually reads off a measured time trace.

44Run the worked example

We call the solver with a=1, b=2, c=5, y(0)=1, y'(0)=0. The characteristic equation is r² + 2r + 5 = 0, giving r = −1 ± 2i. So α = −1 (decay), β = 2 (angular frequency). Plugging the initial conditions gives C₁ = 1, C₂ = 0.5, and A = √1.25 ≈ 1.118.

50Print a sanity check

We print y at t = π/4. By hand: e^(−π/4)(cos(π/2) + 0.5·sin(π/2)) = 0.4559 · (0 + 0.5) ≈ 0.2280. NumPy returns 0.2280 — confirming the code matches the derivation.

38 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_complex_roots(a, b, c, y0, yp0, t):
5    """
6    Solve  a*y'' + b*y' + c*y = 0  in the COMPLEX-ROOTS case
7    (i.e. when the discriminant b^2 - 4ac is negative).
8
9    Returns y(t) and a dict of the physical quantities.
10    """
11    disc = b**2 - 4*a*c
12    assert disc < 0, "complex roots only — got disc >= 0"
13
14    # Roots:  r = alpha ± i beta
15    alpha = -b / (2*a)
16    beta  = np.sqrt(-disc) / (2*a)
17
18    # Apply initial conditions  y(0)=y0,  y'(0)=yp0
19    #   y(0)  = C1                 -> C1 = y0
20    #   y'(0) = alpha*C1 + beta*C2 -> C2 = (yp0 - alpha*y0)/beta
21    C1 = y0
22    C2 = (yp0 - alpha*y0) / beta
23
24    # Real-valued solution
25    y = np.exp(alpha*t) * (C1*np.cos(beta*t) + C2*np.sin(beta*t))
26
27    # Amplitude-phase form  y(t) = A * e^(alpha t) * cos(beta t - phi)
28    A   = np.sqrt(C1**2 + C2**2)
29    phi = np.arctan2(C2, C1)
30
31    return y, {
32        "alpha": alpha,         # decay rate  (real part of the root)
33        "beta":  beta,          # angular freq (imag part of the root)
34        "C1": C1, "C2": C2,
35        "amplitude": A,
36        "phase": phi,
37        "period": 2*np.pi / beta,
38        "half_life": np.log(2)/-alpha if alpha < 0 else float("inf"),
39    }
40
41
42# ---- worked example: y'' + 2 y' + 5 y = 0,  y(0)=1,  y'(0)=0 ----
43t = np.linspace(0, 6, 400)
44y, info = solve_complex_roots(a=1, b=2, c=5, y0=1.0, yp0=0.0, t=t)
45
46print("alpha    =", info["alpha"])
47print("beta     =", info["beta"])
48print("C1, C2   =", info["C1"], info["C2"])
49print("A, phi   =", info["amplitude"], info["phase"])
50print("period   =", info["period"])
51print("y(pi/4)  =", y[int(400*(np.pi/4)/6)])

Why hand-derive when scipy.integrate.solve_ivp exists?

Closed-form solutions let you reason about a system before you simulate it. You can read off the period, the half-life, the steady-state behaviour, all without running a single step. A numerical solver tells you what happens; the closed form tells you why.

PyTorch: The Matrix-Exponential Picture

There is one routine that solves ALL three cases (real distinct, repeated, complex) without any branching on the discriminant: the matrix exponential. Convert the 2nd-order ODE into a first-order linear system $Y' = A\,Y$ with $Y = (y, y')^T$ , and the solution is just $Y(t) = e^{At}\,Y(0)$ .

The eigenvalues of $A$ are exactly the roots of the characteristic equation. For our worked example PyTorch will print $(-1 + 2i,\ -1 - 2i)$ — the same pair we derived by hand.

matrix-exponential solver via torch.linalg.matrix_exp

🐍python

Explanation(8)

Code(39)

1Import PyTorch

We use torch.linalg.matrix_exp, the matrix exponential. For a linear constant-coefficient ODE Y' = A·Y, the solution is literally Y(t) = e^(At)·Y(0) — the matrix exponential IS the propagator. PyTorch ships this as a single call.

11Build the system matrix A

Top row [0, 1] encodes y' = v. Bottom row [-c/a, -b/a] encodes v' = -(c/a)y - (b/a)v. Stacking these turns the 2nd-order equation into a 2-D first-order linear system Y' = A·Y with Y = [y, v].

13Stack the initial conditions

Y(0) = [y(0), y'(0)] — exactly the two initial conditions the 2nd-order ODE requires, now living in a single 2-vector.

16Eigenvalues == characteristic roots

Crucial connection: the eigenvalues of A are exactly the roots r of the characteristic equation a·r² + b·r + c = 0. For our example a=1, b=2, c=5, NumPy/PyTorch prints (-1+2i) and (-1-2i) — the complex conjugate pair we derived by hand. This is why eigenvalue language and root language are interchangeable for linear ODEs.

20Compute the propagator for every t

For each time t_i we form A·t_i (a 2×2 matrix), take its matrix exponential (also 2×2), and multiply by Y(0). The result is a 2-vector [y(t_i), y'(t_i)]. We stack all of them into a (T, 2) tensor.

24Take the y component

Slice column 0 to extract y(t) and drop the velocity. (For phase-plane plots you'd keep both columns.)

30Cross-check against the closed form

We compute y(t) = e^(-t)(cos(2t) + 0.5·sin(2t)) directly and compare. The matrix-exponential method should agree to ~1e-14 in float64 — the two routes give literally the same number. This is the strongest possible sanity check.

34Why this matters

matrix_exp is a single API that solves all three cases (real distinct, repeated, complex) uniformly — no if/else on the discriminant, no special-case formulas. It is also exactly the routine inside torchdiffeq and the linear flow of every Neural ODE. Understanding the 2×2 case scales to N×N coupled oscillators (modal analysis) line-for-line.

31 lines without explanation

1import torch
2
3def solve_via_matrix_exp(a, b, c, y0, yp0, t):
4    """
5    Solve  a*y'' + b*y' + c*y = 0  exactly using the MATRIX exponential.
6
7    Convert to first-order system Y' = A Y with Y = [y, y'].
8    Then  Y(t) = exp(A * t) @ Y(0)  for any t.
9
10    This works for ALL three cases (distinct, repeated, complex) without
11    branching on the discriminant.
12    """
13    A = torch.tensor([[0.0,         1.0],
14                      [-c/a,        -b/a]], dtype=torch.float64)
15    Y0 = torch.tensor([y0, yp0], dtype=torch.float64)
16
17    # eigenvalues of A == roots of the characteristic equation
18    eigvals = torch.linalg.eigvals(A)
19    print("eigenvalues of A =", eigvals)
20
21    # Vectorised: matrix_exp(A * t_i) for every t_i in the grid
22    Y = torch.stack([
23        torch.linalg.matrix_exp(A * ti) @ Y0
24        for ti in t
25    ])
26    return Y[:, 0]   # take only the y component, drop y'
27
28
29# Same worked example: y'' + 2 y' + 5 y = 0
30t = torch.linspace(0.0, 6.0, 400, dtype=torch.float64)
31y_pt = solve_via_matrix_exp(1.0, 2.0, 5.0, 1.0, 0.0, t)
32
33# Hand-derived closed form for cross-check
34alpha, beta = -1.0, 2.0
35C1, C2 = 1.0, 0.5
36y_exact = torch.exp(alpha*t) * (C1*torch.cos(beta*t) + C2*torch.sin(beta*t))
37
38err = torch.max(torch.abs(y_pt - y_exact)).item()
39print("max abs error vs closed-form:", err)

Eigenvalues = roots: the unifying view

The characteristic equation $a r^2 + b r + c = 0$ is, up to sign, exactly the characteristic polynomial of the 2×2 system matrix $A = \begin{pmatrix} 0 & 1 \\ -c/a & -b/a \end{pmatrix}$ . So "find the roots of the characteristic equation" and "find the eigenvalues of $A$ " are the same computation. This connection scales perfectly: an $n$ -th order ODE becomes an $n \times n$ first-order system, and the eigenvalues still classify the dynamics.

Where Complex Roots Show Up in the Real World

🪨 Car suspension

A shock absorber over a bump: $m\ddot{x} + c\dot{x} + kx = 0$ . Designers pick $c$ so that $\zeta \approx 0.7$ — fast settle with only a barely-perceptible overshoot.

⚡ RLC circuit

$L\ddot{q} + R\dot{q} + q/C = 0$ . Same equation. $R$ plays the role of damping; a low-R resonator (e.g. a quartz crystal) has tiny $\alpha$ and rings for millions of cycles.

🧬 Molecular vibrations

A diatomic molecule's stretching mode is a damped oscillator whose $\beta$ tells you the infrared-absorption frequency. Spectroscopy = measuring $\beta$ .

🤖 Gradient descent with momentum

Heavy-ball / momentum SGD near a quadratic minimum behaves exactly like an underdamped oscillator. Picking the momentum coefficient badly puts you in the underdamped regime — you see the loss oscillate. Pick it as if you were tuning $\zeta \approx 1$ for critical damping.

🏙️ Building sway in wind

Skyscrapers have tuned mass dampers (huge pendulums in the penthouse) sized to push the building's $\zeta$ up to a comfortable level. Without them, residents feel underdamped sway for tens of seconds after a gust.

🎛️ Audio synthesizers

A plucked-string sound is a damped oscillator chosen to give a long ringing tail. The decay rate $|\alpha|$ sets "sustain"; $\beta$ sets pitch.

Common Pitfalls

1. Confusing β with f (Hz)

$\beta$ is in radians per second. The ordinary frequency in Hz is $f = \beta/(2\pi)$ . Period is $T = 2\pi/\beta = 1/f$ . Mixing these up by a factor of $2\pi$ is the most common numerical bug in textbook problems.

2. Forgetting that β is √(−Δ)/(2a), not √Δ/(2a)

In the complex branch $\Delta < 0$ , so $\sqrt{\Delta}$ is not real. The correct imaginary part uses $\sqrt{-\Delta}$ . If you blindly type sqrt(disc) into Python and it returns a NaN, you have just rediscovered this trap.

3. Treating the envelope as the solution

$A\,e^{\alpha t}$ bounds the solution but is not itself a solution. The actual $y(t)$ touches the envelope only at the peaks of the cosine. Average behaviour: zero crossings every $\pi/\beta$ .

4. atan vs atan2 for the phase

$\arctan(C_2/C_1)$ alone can't tell $(C_1, C_2) = (1, 1)$ from $(-1, -1)$ — both give ratio 1. Always use $\mathrm{atan2}(C_2, C_1)$ when computing phase, exactly as the Python code above does.

5. Numerical issues near ζ = 1

Right at the critical damping boundary the formulas split. Code that branches on $\Delta < 0$ vs $\Delta = 0$ needs a tolerance — pure equality almost never triggers. The matrix-exponential approach sidesteps this completely since it has no branches.

Summary

What to take away

Complex roots happen iff $\Delta = b^2 - 4ac < 0$ .
The roots are $r = \alpha \pm i\beta$ with $\alpha = -b/(2a)$ and $\beta = \sqrt{-\Delta}/(2a)$ .
Euler's formula $e^{i\theta} = \cos\theta + i\sin\theta$ turns the complex exponential into a sum of sines and cosines.
The real general solution is $y(t) = e^{\alpha t}\bigl(C_1\cos\beta t + C_2\sin\beta t\bigr)$ , equivalently $y(t) = A\,e^{\alpha t}\cos(\beta t - \varphi)$ .
$\alpha$ is the decay rate of the envelope; $\beta$ is the angular frequency of the oscillation. Left half-plane = stable; imaginary axis = perpetual motion; right half-plane = unstable growth.
Engineers reparameterise as $(\omega_n, \zeta)$ . The damping ratio $\zeta$ classifies all four regimes: undamped (0), underdamped (0–1), critical (1), overdamped (>1).
The same equation is the textbook model for springs, RLC circuits, molecules, momentum-SGD, and skyscrapers in the wind. Master this section and you have mastered every linear oscillator in physics and engineering.

Looking ahead: next section we close the loop with the third case — the repeated root — and explain why the second independent solution suddenly grows a factor of $t$ . After that we leave the homogeneous world entirely and add a forcing term to the right-hand side.