Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Recognize a nonhomogeneous second-order linear ODE $\;\; ay'' + by' + cy = g(t)$ and explain why its solution splits as $y = y_h + y_p$ .
Choose the correct trial form for $y_p$ when the forcing $g(t)$ is a polynomial, exponential, sinusoid, or product/sum of these.
Apply the "multiply by $t$ " rule when the trial form overlaps the homogeneous solution (the resonance case).
Determine the unknown coefficients by substituting the trial back into the ODE and matching like terms.
Apply initial conditions correctly to the total $y = y_h + y_p$ , not to $y_h$ alone.
Implement the method in Python and verify it with PyTorch autograd.

The Big Picture: Forcing a System

"Section 22.01 showed how a system moves when nothing pushes it. This section asks: what happens when something does?"

Until now every ODE we have written had a zero on the right-hand side. Physically that meant a system left alone: a swinging pendulum after you let go, a charged capacitor discharging through a resistor, a spring released from its stretched position. The motion was governed entirely by the system's own internal dynamics — the homogeneous equation.

But the most useful problems involve external input: a constant gravitational pull, a sinusoidal voltage source, an oscillating piston, a sudden hammer-strike. Those become a nonhomogeneous ODE:

a\,y'' + b\,y' + c\,y = g(t).

The function $g(t)$ on the right-hand side is the forcing term — the input that drives the system. Solving this ODE answers a single, very practical question:

Given an external push $g(t)$ , what motion does the system produce?

The Core Idea (in one line)

For a linear ODE, the answer always has the form

y(t) = \underbrace{y_h(t)}_{\text{free response}} \;+\; \underbrace{y_p(t)}_{\text{forced response}}.

$y_h$ is what the system would do if nothing pushed it (the homogeneous solution from Section 22.01–03). $y_p$ is one specific motion that the ODE produces in response to the forcing. Adding the two gives every possible solution.

The Superposition Principle

Why can we split the solution into two pieces? Because the operator $L[y] = ay'' + by' + cy$ is linear. Linearity says: for any constants $\alpha, \beta$ and any twice-differentiable $u, v$ ,

L[\alpha u + \beta v] = \alpha\,L[u] + \beta\,L[v].

Now suppose $y_p$ is any single solution of $L[y] = g$ , and let $y_h$ be the general solution of the homogeneous equation $L[y] = 0$ . Then:

L[y_h + y_p] = L[y_h] + L[y_p] = 0 + g = g. \;\checkmark

So $y_h + y_p$ is also a solution. Conversely, if $Y$ is any solution then $Y - y_p$ satisfies $L[Y - y_p] = g - g = 0$ , meaning $Y - y_p$ is a homogeneous solution — call it $y_h$ . So $Y = y_h + y_p$ . Every solution decomposes this way.

The strategy

To solve $L[y] = g$ we only need:

The general homogeneous solution $y_h = C_1 y_1 + C_2 y_2$ (sections 22.01–22.03).
One particular solution $y_p$ — any single function that satisfies the full equation.
Add them and fit $C_1, C_2$ to the initial conditions on the TOTAL $y$ .

Finding $y_p$ is what this section is about.

The Method of Undetermined Coefficients

For certain special forcings $g(t)$ — polynomials, exponentials, sines, cosines, and their sums and products — there is a remarkably clean recipe:

Guess a trial form for $y_p$ with the same algebraic "family" as $g(t)$ , but with unknown coefficients.
Substitute the trial into the ODE.
Match coefficients of like terms on both sides.
Solve the resulting linear system for the unknowns.

Why does this work? Because differentiation maps these families to themselves. The derivative of a polynomial is a polynomial. The derivative of $e^{\alpha t}$ is $\alpha e^{\alpha t}$ — still an exponential of the same rate. The derivative of a sinusoid of frequency $\omega$ is a sinusoid of frequency $\omega$ . So the LHS of the ODE stays inside the family of the RHS, and the matching procedure terminates after finitely many algebraic steps.

What "undetermined" means

The coefficients ( $A, B, \ldots$ ) in the trial form are unknown until step 3 forces them. They are "determined" only AFTER you substitute and match — hence the name. There is no calculus magic; it is bookkeeping.

Catalog of Trial Forms

Memorize this table. It is the whole method, on one page.

Forcing g(t)	Trial form y_p	Why
k (constant)	A	Derivatives of a constant are 0, so plug A into ODE: cA = k.
p₀ + p₁ t (linear)	A + B t	Same degree as the forcing. Polynomial in, polynomial out.
Polynomial of degree n	A₀ + A₁t + … + Aₙ tⁿ	Differentiation lowers degree; the highest-degree term of y_p must match the forcing.
K e^(αt)	A e^(αt)	e^(αt) is an eigenfunction of d/dt — multiplying by a constant is the only operation needed.
K cos(ωt) or K sin(ωt)	A cos(ωt) + B sin(ωt)	You need BOTH; differentiating cos produces sin, only the pair is closed under d/dt.
K e^(αt) cos(ωt)	e^(αt) (A cos(ωt) + B sin(ωt))	Product of an exponential and a sinusoid — combine both rules.
Sum of the above	Corresponding sum of trial forms	Linearity: solve for each piece, then add.

The matching principle

Always write the trial form so it contains every term you could get by differentiating it twice. That is why the sinusoid trial has both $\cos$ and $\sin$ : even if the forcing is just $\cos$ , the LHS of the ODE will produce sine terms.

The Resonance Rule: Multiply by $t$

The catalog above breaks down in one important case: when the trial form is already a solution of the homogeneous equation.

Consider $y'' - 3y' + 2y = e^{t}$ . The homogeneous roots are $r = 1, 2$ , so $y_h = C_1 e^{t} + C_2 e^{2t}$ . Try $y_p = A e^{t}$ :

A e^{t} - 3A e^{t} + 2A e^{t} = 0 \cdot A e^{t} = 0 \neq e^{t}.

The LHS collapses to zero no matter what $A$ we choose — because $e^{t}$ is already in the kernel of $L$ . Algebraically: $Q(1) = 1 - 3 + 2 = 0$ , so the formula $A = K / Q(\alpha)$ divides by zero.

The Resonance Rule

If the trial form overlaps the homogeneous solution $y_h$ , multiply the trial by $t$ . If it still overlaps (when the root is repeated), multiply by $t^2$ .

In symbols: the "multiplicity bump" equals the algebraic multiplicity of the root that causes the clash.

Situation	Bumped trial form
α is NOT a root of the characteristic equation	A e^(αt)
α is a SIMPLE root	A t · e^(αt)
α is a DOUBLE root	A t² · e^(αt)
Forcing K cos(ωt), iω is NOT a homogeneous root	A cos(ωt) + B sin(ωt)
Forcing K cos(ωt), iω IS a homogeneous root (undamped at natural freq)	t · (A cos(ωt) + B sin(ωt))

Where the t comes from

This is exactly the same phenomenon you saw in Section 22.03 with $t\,e^{rt}$ : when a root coalesces, the second linearly independent solution always gains a factor of $t$ . Resonance is the SAME algebra — a root of the characteristic equation has "moved into" the forcing, so the same fix applies.

Physically, resonance is the most striking phenomenon in linear systems: a sinusoidal force at the system's natural frequency builds the amplitude without bound. This is why opera singers shatter wineglasses, why platoons break step on bridges, and why every engineer designs damping into rotating machinery.

Worked Example (Step-by-Step)

Solve $y'' - 5y' + 6y = 4 e^{4t}$ subject to $y(0) = 0,\; y'(0) = 0$ .

Click to expand the full hand calculation

Step 1 — Solve the homogeneous equation first.

Char. eq: $r^2 - 5r + 6 = 0$ . Factor: $(r-2)(r-3) = 0$ so $r_1 = 2,\, r_2 = 3$ . Hence

y_h = C_1 e^{2t} + C_2 e^{3t}.

Step 2 — Check the forcing for resonance.

Forcing is $g(t) = 4 e^{4t}$ , an exponential with rate $\alpha = 4$ . Is $4$ a root of the characteristic equation? $Q(4) = 16 - 20 + 6 = 2 \neq 0$ . No clash — use the standard trial form.

Step 3 — Pick the trial form.

y_p = A\,e^{4t}.

Step 4 — Differentiate.

$y_p' = 4A e^{4t},\quad y_p'' = 16 A e^{4t}.$

Step 5 — Substitute into the ODE.

16 A e^{4t} - 5 \cdot 4 A e^{4t} + 6 A e^{4t} = 4 e^{4t}.

Combine: $(16 - 20 + 6) A e^{4t} = 2 A e^{4t} = 4 e^{4t}$ .

Step 6 — Solve for A.

$2A = 4 \Rightarrow A = 2$ . So $y_p = 2 e^{4t}$ .

Step 7 — Write the general solution.

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

Step 8 — Apply $y(0) = 0$ .

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

Step 9 — Apply $y'(0) = 0$ .

$y'(t) = 2 C_1 e^{2t} + 3 C_2 e^{3t} + 8 e^{4t}$ so $y'(0) = 2 C_1 + 3 C_2 + 8 = 0$ $\Rightarrow$ $2 C_1 + 3 C_2 = -8$ .

Step 10 — Solve the 2×2 system.

From the first equation, $C_1 = -2 - C_2$ . Substitute: $2(-2 - C_2) + 3 C_2 = -8$ $\Rightarrow$ $-4 + C_2 = -8$ $\Rightarrow$ $C_2 = -4$ . Then $C_1 = -2 - (-4) = 2$ .

Step 11 — Particular solution.

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

Step 12 — Sanity check.

At $t=0$ : $y(0) = 2 - 4 + 2 = 0$ ✓. $y'(0) = 4 - 12 + 8 = 0$ ✓.

At $t = 0.5$ : $y(0.5) = 2 e^1 - 4 e^{1.5} + 2 e^2 \approx 5.437 - 17.925 + 14.778 \approx 2.290$ . The Python solver below returns $y(0.5) \approx 2.288$ , matching to 3 decimals.

Step 13 — Notice the cancellation.

The particular piece alone is $y_p(0.5) = 2 e^2 \approx 14.78$ , but the TOTAL $y(0.5) \approx 2.29$ — much smaller. The homogeneous piece is doing real work, fighting the forcing in the early transient so the IC's $y(0) = y'(0) = 0$ are honoured. Many students confuse $y_p$ with the actual physical motion — it is not, by itself.

Interactive Forcing Explorer

Pick a forcing type — polynomial, exponential, or sinusoid — and watch four curves update simultaneously: the forcing $g(t)$ itself (dashed gray), the transient $y_h$ (green), the particular solution $y_p$ (blue), and the total $y = y_h + y_p$ (red). Try every preset.

Loading interactive explorer…

Three things to discover

Switch to exponential forcing and slide $\alpha$ through the homogeneous roots (with the default $a=1, b=3, c=2$ they are $-1$ and $-2$ ). At those exact points an amber warning appears: the trial gets multiplied by $t$ .
Switch to sinusoid, set $a=1, b=0, c=4$ , and sweep $\omega$ . The natural frequency is $\omega_0 = 2$ . Park $\omega$ exactly at 2 — the response grows like $t \sin(2t)$ , the textbook resonance.
Add even tiny damping $b = 0.3$ and re-do the sweep. The unbounded growth disappears; the steady amplitude becomes large but finite. This is exactly why suspension bridges have dampers.

Transient vs Steady State

When the homogeneous solution decays — $\mathrm{Re}(r_1), \mathrm{Re}(r_2) < 0$ , the stable case — the picture splits naturally into two regimes.

Regime	Dominated by	Time scale
Transient	y_h (initial conditions, free response)	Decays on the time-scale 1 / \|Re(r)\|.
Steady state	y_p (the forcing response)	Lasts forever — same shape as g(t).

After the transient dies, $y(t) \approx y_p(t)$ . Engineers call this the steady-state response and spend enormous effort designing it: amplifier frequency response, loudspeaker damping, automotive suspension. The transient is the "ringing" before the system settles.

Why the initial conditions usually do not matter long-term

Notice that the IC's only show up in $y_h$ (they pin $C_1, C_2$ ). If $y_h \to 0$ , the IC information washes out. A stable, forced linear system forgets its starting conditions. This is the rigorous reason a hi-fi system doesn't depend on how you switched it on.

Computation in Python

We turn the recipe into code. Given $(a, b, c)$ and a description of $g(t)$ , the function below returns the total $y$ , the transient $y_h$ , and the steady-state piece $y_p$ on any time grid. It dispatches on the forcing type, builds the trial-form coefficients in closed form, applies the resonance bump where needed, then solves a 2×2 system for the homogeneous constants from the total-y initial conditions.

Analytic solver for three families of forcings

Trial-form dispatch + resonance handling + IC fitting

🐍undetermined_coeffs.py

Explanation(15)

Code(116)

1Import NumPy

NumPy lets us evaluate the closed-form solution on an entire time grid in one expression. No Python-level loop is needed.

3Function header

We take the ODE coefficients a, b, c, a description of the forcing g(t), the two initial conditions, and a time grid. A second-order ODE needs TWO initial conditions because the homogeneous solution carries two free constants C₁, C₂.

23Read the forcing type

The dispatcher pattern: we switch on forcing['kind']. Each kind has its own trial form and its own algebra for the undetermined coefficient. The user does not have to know the formulas — they hand us a dict and we look up the right rule.

25Exponential branch

For g(t) = K e^(αt) the trial is y_p = A e^(αt). Substituting: a·A α² + b·A α + c·A = K. The bracket is the characteristic polynomial Q(α). So A = K / Q(α). The whole method collapses to one division.

27Check the resonance condition

If Q(α) is essentially zero, α IS a root of the characteristic equation — the trial e^(αt) lives inside the homogeneous solution. Dividing by zero would blow up, and physically we are forcing the system at its natural rate. The fix below: multiply the trial by t.

33Resonant exponential trial

Substituting y_p = A·t·e^(αt) and using y_p'' = A(α²t + 2α)e^(αt) gives A·(2aα + b) e^(αt) = K e^(αt). The bracket 2aα + b is Q′(α) — the derivative of the characteristic polynomial at α. For a SIMPLE root Q′(α) ≠ 0 so A is well defined.

37Sinusoidal branch

For g(t) = K cos(ωt) the natural trial includes BOTH cosine and sine: y_p = A cos(ωt) + B sin(ωt). You always need both because differentiating a cosine produces a sine — only the combination is closed under the ODE.

39Build the 2×2 linear system

Matching coefficients of cos and sin gives: [D, E; -E, D] · [A; B] = [K; 0], where D = c − aω² (the 'stiffness mismatch') and E = bω (the 'damping torque'). Determinant D² + E² is non-negative.

41Solve for A and B

By Cramer's rule, A = K·D / (D² + E²) and B = K·E / (D² + E²). The amplitude √(A² + B²) = K / √(D² + E²) is the classic resonance curve: it peaks where D ≈ 0 (i.e. ω ≈ √(c/a)) and is suppressed by damping E = bω.

47Pure-resonance branch (b = 0, ω² = c/a)

When the system is undamped AND forced exactly at its natural frequency, both D and E are zero. The trial must be y_p = t·(A cos ωt + B sin ωt). Substitution gives A = 0 and B = K/(2aω). The solution grows linearly in t — unbounded resonance.

56Polynomial branch

For g(t) = p₀ + p₁t the trial is the SAME degree polynomial: y_p = A + B·t. Substituting (note y_p'' = 0) gives c·B = p₁ and b·B + c·A = p₀. The if-branch handles the usual case c ≠ 0. The c = 0 case is in the page text (multiply by t).

70Discriminant for the homogeneous part

We need y_h = C₁ y₁ + C₂ y₂. The three branches (distinct real / complex / repeated) come from the sign of Δ = b² − 4ac — exactly the same classifier we built in Section 22.01. The free constants C₁, C₂ are pinned later by the initial conditions on the TOTAL y.

79Initial conditions on the homogeneous part

Critical step: the IC's y(0)=y0, y'(0)=yp0 apply to the TOTAL y = y_h + y_p. So y_h(0) = y0 − y_p(0) and y_h'(0) = yp0 − y_p'(0). Many students apply the IC's directly to y_h and get the wrong answer.

99Sum on the grid

Once both pieces are vectorized over t, we add them elementwise. The transient yh decays away (if Re(roots) < 0), leaving y ≈ yp at large t — the steady state.

106Run the worked example

y'' − 5y' + 6y = 4e^(4t) with y(0)=y'(0)=0. Homogeneous roots are 2, 3; 4 is NOT a root, so the trial is A e^(4t) and A = 4/(16 − 20 + 6) = 2. Total: y = 2e^(2t) − 4e^(3t) + 2e^(4t). At t=0.5 the particular piece is 2·e^2 ≈ 14.778; the total is much smaller (2.288) because the homogeneous part nearly cancels it.

101 lines without explanation

1import numpy as np
2
3def solve_nonhomogeneous(a, b, c, forcing, y0, yp0, t):
4    """
5    Solve  a*y'' + b*y' + c*y = g(t)  by:
6      1) finding a particular solution y_p with undetermined coefficients
7      2) finding y_h (general homogeneous solution)
8      3) fitting C1, C2 of y_h so y(0)=y0 and y'(0)=yp0
9
10    Parameters
11    ----------
12    a, b, c  : float       constant coefficients
13    forcing  : dict        {"kind": "exp", "K": float, "alpha": float}
14                    or     {"kind": "cos", "K": float, "omega": float}
15                    or     {"kind": "poly", "p0": float, "p1": float}
16    y0, yp0  : float       initial conditions for the TOTAL y
17    t        : ndarray     time grid
18
19    Returns
20    -------
21    y_total, y_h, y_p   each a 1-D numpy array on the grid t.
22    """
23
24    # ---------- 1. Build the particular solution ----------
25    kind = forcing["kind"]
26
27    if kind == "exp":
28        K, alpha = forcing["K"], forcing["alpha"]
29        Q = a * alpha**2 + b * alpha + c       # char polynomial at alpha
30        if abs(Q) > 1e-9:
31            A = K / Q
32            y_p  = lambda u: A * np.exp(alpha * u)
33            y_pd = lambda u: A * alpha * np.exp(alpha * u)
34        else:
35            # alpha IS a root  →  bump by t
36            Qp = 2 * a * alpha + b
37            A = K / Qp
38            y_p  = lambda u: A * u * np.exp(alpha * u)
39            y_pd = lambda u: A * (1 + alpha * u) * np.exp(alpha * u)
40
41    elif kind == "cos":
42        K, w = forcing["K"], forcing["omega"]
43        D, E = c - a * w**2, b * w              # 2x2 linear system
44        det = D**2 + E**2
45        if det > 1e-9:
46            A = K * D / det
47            B = K * E / det
48            y_p  = lambda u: A * np.cos(w * u) + B * np.sin(w * u)
49            y_pd = lambda u: -A * w * np.sin(w * u) + B * w * np.cos(w * u)
50        else:
51            # Pure resonance:  b=0 and w^2 = c/a
52            A, B = 0.0, K / (2 * a * w)
53            y_p  = lambda u: A * u * np.cos(w * u) + B * u * np.sin(w * u)
54            y_pd = lambda u: (A + B * w * u) * np.cos(w * u) \
55                           + (B - A * w * u) * np.sin(w * u)
56
57    elif kind == "poly":
58        p0, p1 = forcing["p0"], forcing["p1"]
59        if abs(c) > 1e-9:
60            B = p1 / c
61            A = (p0 - b * B) / c
62            y_p  = lambda u: A + B * u
63            y_pd = lambda u: B * np.ones_like(u)
64        else:
65            raise NotImplementedError("c=0 polynomial case handled in the page.")
66
67    else:
68        raise ValueError(f"unknown forcing kind: {kind}")
69
70    # ---------- 2. Homogeneous part: solve characteristic eq. ----------
71    disc = b**2 - 4 * a * c
72    if disc > 1e-9:
73        sqrtD = np.sqrt(disc)
74        r1 = (-b + sqrtD) / (2 * a)
75        r2 = (-b - sqrtD) / (2 * a)
76        y_h_basis = lambda u, C1, C2: C1 * np.exp(r1 * u) + C2 * np.exp(r2 * u)
77        # y_h'(t) = C1 r1 e^{r1 t} + C2 r2 e^{r2 t}
78        # IC for y_h: y_h(0) = y0 - y_p(0),  y_h'(0) = yp0 - y_p'(0)
79        Y0  = y0  - y_p(np.array(0.0))
80        YP0 = yp0 - y_pd(np.array(0.0))
81        C2 = (YP0 - r1 * Y0) / (r2 - r1)
82        C1 = Y0 - C2
83        y_h_eval = lambda u: C1 * np.exp(r1 * u) + C2 * np.exp(r2 * u)
84
85    elif disc < -1e-9:
86        alpha_h = -b / (2 * a)
87        beta_h  = np.sqrt(-disc) / (2 * a)
88        Y0  = y0  - y_p(np.array(0.0))
89        YP0 = yp0 - y_pd(np.array(0.0))
90        C1 = Y0
91        C2 = (YP0 - alpha_h * Y0) / beta_h
92        y_h_eval = lambda u: np.exp(alpha_h * u) * (C1 * np.cos(beta_h * u)
93                                                    + C2 * np.sin(beta_h * u))
94    else:
95        r = -b / (2 * a)
96        Y0  = y0  - y_p(np.array(0.0))
97        YP0 = yp0 - y_pd(np.array(0.0))
98        C1 = Y0
99        C2 = YP0 - r * Y0
100        y_h_eval = lambda u: (C1 + C2 * u) * np.exp(r * u)
101
102    # ---------- 3. Sum on the time grid ----------
103    yh = y_h_eval(t)
104    yp = y_p(t)
105    return yh + yp, yh, yp
106
107
108# -------- worked example: y'' - 5y' + 6y = 4 e^{4t}, y(0)=0, y'(0)=0 --------
109t = np.linspace(0, 1.2, 200)
110y, yh, yp = solve_nonhomogeneous(
111    a=1, b=-5, c=6,
112    forcing={"kind": "exp", "K": 4, "alpha": 4},
113    y0=0, yp0=0, t=t,
114)
115print("y_p(0.5) =", yp[83])       # should match 2 * exp(2) ≈ 14.778
116print("y(0.5)   =", y[83])

The resonance amplitude curve

For sinusoidal forcing the steady-state amplitude as a function of $\omega$ is one of the most studied curves in engineering — every audio EQ, every seismometer's response, every NMR linewidth comes from it. Plot it once and remember the shape.

Steady-state amplitude vs forcing frequency

🐍resonance_curve.py

Explanation(6)

Code(28)

1Import NumPy

Resonance curves are scalar-in/scalar-out; vectorizing over the ω grid lets us plot the whole curve with no loop.

3Function: amplitude as a function of ω

We computed A and B for the cosine trial in the main solver. The amplitude of the steady-state oscillation is √(A² + B²) = K / √(D² + E²). This single formula contains everything: the height of the resonance peak, its location, and how damping flattens it.

10Build D and E vectorized

D = c − aω² is the 'stiffness mismatch': zero exactly at the undamped natural frequency ω₀ = √(c/a). E = bω is the 'damping reaction' that grows linearly with ω. Both are arrays the same shape as omega_grid.

14Set parameters and the ω grid

a=1, c=4 gives ω₀ = 2. Tiny damping b=0.3 produces a moderate peak; b=0 would produce an infinite peak (true resonance) and b=2 would smear the peak away. K=1 is a unit forcing.

19Find the resonance peak numerically

np.argmax returns the index of the largest amplitude. With small damping it sits just below ω₀. Theory: the peak is at √(c/a − (b/(2a))²) = √(4 − 0.0225) ≈ 1.994.

24Print the diagnostics

Three frequencies to keep straight: (i) ω₀ = √(c/a) — undamped natural; (ii) ω_d = √(ω₀² − (b/2a)²) — damped natural; (iii) the peak of the amplitude response — slightly different from ω_d (it is √(ω₀² − 2(b/2a)²) for amplitude). Engineers use them interchangeably when damping is light.

22 lines without explanation

1import numpy as np
2
3def amplitude_response(a, b, c, K, omega_grid):
4    """
5    For a sinusoidal forcing K cos(omega t), the steady-state amplitude is
6        |y_p|_max = K / sqrt(D^2 + E^2),     D = c - a w^2,  E = b w.
7
8    This is the famous resonance curve. Returns one value per omega.
9    """
10    D = c - a * omega_grid**2
11    E = b * omega_grid
12    return K / np.sqrt(D**2 + E**2)
13
14a, b, c, K = 1.0, 0.3, 4.0, 1.0
15omegas = np.linspace(0.1, 4.0, 400)
16amps   = amplitude_response(a, b, c, K, omegas)
17
18# Find the omega that maximizes amplitude
19i_star  = np.argmax(amps)
20w_peak  = omegas[i_star]
21w_nat   = np.sqrt(c / a)          # undamped natural frequency
22w_damped = np.sqrt(c / a - (b / (2 * a))**2)   # damped natural frequency
23
24print(f"natural omega_0 = {w_nat:.3f}")
25print(f"damped resonance peak omega ~ {w_peak:.3f}")
26print(f"theory damped omega = {w_damped:.3f}")
27print(f"peak amplitude       = {amps[i_star]:.3f}")
28print(f"DC amplitude (w=0)   = {amps[0]:.3f}")

PyTorch View: Autograd Verification

Hand-derived ODE solutions are easy to get wrong. The cheapest test we know: feed the candidate back into the ODE and check the residual $R(y) = a y'' + b y' + c y - g(t)$ is machine-zero. PyTorch's autograd lets us compute $y'$ and $y''$ automatically — no chain-rule by hand, no finite differences.

Verify the closed form by automatic differentiation

🐍residual_autograd.py

Explanation(11)

Code(39)

1Import PyTorch

PyTorch's autograd lets us differentiate any tensor-valued Python function symbolically (well, automatically) instead of by hand. That makes ODE residual checking almost free.

3Higher-order factory: build a residual

We want to check 'does my candidate y(t) satisfy the ODE?' for any forcing g and any closed form. The factory captures a, b, c, g in a closure and returns a function residual(t, y_fn) that produces R = a·y'' + b·y' + c·y − g(t).

8Detach + require gradients

We rebuild t as a leaf tensor with requires_grad=True so that downstream computations track t in the autograd graph. Without this, torch.autograd.grad would have nothing to differentiate with respect to.

9Evaluate the candidate y(t)

y_fn is any differentiable function. Here we will pass in the closed-form 2·e^(2t) − 4·e^(3t) + 2·e^(4t), but it could equally be a neural network or a polynomial.

10First derivative via autograd

torch.autograd.grad(y.sum(), t) returns ∂(Σy)/∂t. Because y depends on t elementwise, this equals ∂y_i/∂t_i = y′(t_i) for each grid point — exactly what we want. create_graph=True keeps the graph alive so we can take a SECOND derivative.

11Second derivative

Same trick again on yd. The result is y''(t) on the grid. The chain create_graph=True → another autograd.grad call is how arbitrary-order derivatives are computed in PyTorch.

14Build the residual

Plug everything into a·y'' + b·y' + c·y − g(t). If y_fn really is a solution, this should be zero (up to float64 round-off). If it is not, the residual tells you where and by how much you missed.

19Define the closed-form candidate

From the worked example: homogeneous part 2e^(2t) − 4e^(3t) (fitted to IC's) + particular 2e^(4t) (from undetermined coefficients). One torch.exp call per term so it is autodifferentiable.

22Define the forcing

g(t) = 4 e^(4t). PyTorch evaluates this as a vector when given a vector t.

27Compute the residual on a fine grid

Run the residual on 200 points in [0, 1.2]. The max absolute residual should be ~1e-13 in float64 — essentially machine epsilon times the magnitudes involved. That is the strongest possible end-to-end sanity check of the derivation.

29Read off y(0) and y′(0)

Confirms both initial conditions are satisfied: y(0) = 2 − 4 + 2 = 0 ✓ and y′(0) = 4 − 12 + 8 = 0 ✓. If either is wrong, you mis-applied the IC's to the WHOLE y (a very common error — see the pitfalls).

28 lines without explanation

1import torch
2
3def make_residual(a, b, c, forcing):
4    """
5    Build the ODE residual  R(y) = a*y'' + b*y' + c*y - g(t),
6    using PyTorch autograd to compute y' and y'' from any differentiable y.
7    """
8    def residual(t, y_fn):
9        t = t.detach().clone().requires_grad_(True)
10        y   = y_fn(t)
11        yd  = torch.autograd.grad(y.sum(),  t, create_graph=True)[0]
12        ydd = torch.autograd.grad(yd.sum(), t, create_graph=True)[0]
13        g = forcing(t)
14        return a * ydd + b * yd + c * y - g
15    return residual
16
17
18# Worked example:  y'' - 5y' + 6y = 4 e^{4t},  y(0)=0, y'(0)=0.
19# Closed form: y(t) = 2 e^{2t} - 4 e^{3t} + 2 e^{4t}.
20
21def y_closed(t):
22    return 2 * torch.exp(2*t) - 4 * torch.exp(3*t) + 2 * torch.exp(4*t)
23
24def forcing(t):
25    return 4 * torch.exp(4 * t)
26
27R = make_residual(a=1.0, b=-5.0, c=6.0, forcing=forcing)
28
29t = torch.linspace(0.0, 1.2, 200, dtype=torch.float64)
30r = R(t, y_closed)
31print("max |residual|        =", r.abs().max().item())
32print("y_closed(0)           =", y_closed(torch.tensor(0.0)).item())
33print("d/dt y_closed at t=0  =", (
34    torch.autograd.grad(
35        y_closed(t.detach().clone().requires_grad_(True)).sum(),
36        t.detach().clone().requires_grad_(True),
37        create_graph=False
38    )[0][0]
39).item())

Connection to physics-informed neural networks (PINNs)

Exactly this residual is what PINNs minimize: they parametrise $y(t; \theta)$ as a neural network and train $\theta$ to make $R(y)$ small at sampled collocation points. The method of undetermined coefficients gives the network the "target" — for the special class of linear ODEs with structured forcings, we already know the closed form, so PINNs should be able to recover it. This section's formulas are the analytic baseline against which every learned solver is measured.

Common Pitfalls

Applying the IC's to y_h instead of y

The initial conditions $y(0) = y_0$ and $y'(0) = y_0'$ apply to the total $y = y_h + y_p$ , not to $y_h$ alone. The correct IC's for $y_h$ are $y_h(0) = y_0 - y_p(0)$ and $y_h'(0) = y_0' - y_p'(0)$ . Forgetting to subtract $y_p(0)$ is the single most common error in this section.

Forgetting the resonance bump

Always check whether your trial form is already in the homogeneous solution. If you skip this check and divide by the characteristic polynomial at $\alpha$ , you will divide by zero (or worse, get a finite but completely wrong answer when round-off masks the zero).

Cos-only trial for cos-only forcing

Even if $g(t) = K \cos(\omega t)$ contains no sine, the trial form MUST include both: $y_p = A\cos\omega t + B\sin\omega t$ . Differentiation will produce sine terms whether you like it or not. A pure-cosine trial leads to an over-determined system with no solution.

Trying undetermined coefficients on g(t) = 1/t or g(t) = tan t

The method only works for the families listed in the catalog: polynomial × exponential × sinusoid. For $1/t, \tan t, \ln t, \sqrt t$ , etc., the space of derivatives never closes — there is no finite trial form. For those forcings, use Variation of Parameters (Section 22.05).

Mistaking y_p for the physical motion

$y_p$ is ONE solution, often without the right initial values. The actual motion is $y_h + y_p$ with $C_1, C_2$ fitted to the IC's. The worked example shows this dramatically: $y_p(0.5) \approx 14.78$ but the actual $y(0.5) \approx 2.29$ .

Summary

Solving $a y'' + b y' + c y = g(t)$ by undetermined coefficients is a four-step recipe:

Solve the homogeneous equation. Get $y_h = C_1 y_1 + C_2 y_2$ and the characteristic roots $r_1, r_2$ .
Choose a trial form for $y_p$ matching the family of $g(t)$ . If the trial overlaps $y_h$ , multiply by $t$ (or $t^2$ if the root is repeated).
Substitute the trial into the ODE and match coefficients of like terms. Solve for the unknown constants.
Write $y = y_h + y_p$ and apply the initial conditions to the TOTAL $y$ to pin $C_1, C_2$ .

Trial-form cheat sheet

Forcing g(t)	Standard trial	Bump if resonant
Polynomial deg n	A₀ + A₁t + … + Aₙtⁿ	Multiply by t (or t²) if r=0 is a root
K e^(αt)	A e^(αt)	Multiply by t (or t²) if α is a root
K cos(ωt) or K sin(ωt)	A cos(ωt) + B sin(ωt)	Multiply by t if ±iω is a root
e^(αt)·polynomial	e^(αt) · (A₀ + … + Aₙtⁿ)	Multiply by t^k where k = multiplicity of α

The slogan

"Guess the same family as the forcing — and if the guess is already in y_h, multiply by t."

Coming next: Section 22.05 introduces Variation of Parameters — a method that handles any continuous forcing $g(t)$ , including $\tan t$ , $\sec t$ , and arbitrary tabulated data. The price you pay is two integrals. Undetermined coefficients is faster when it applies; variation of parameters always applies.