Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Derive the series-RLC differential equation $L\,q'' + R\,q' + \tfrac{1}{C}\,q = 0$ directly from Kirchhoff's voltage law.
Recognize that the RLC equation is the same second-order linear ODE as the mass–spring–damper, so everything you learned in Section 22.06 applies here verbatim.
Classify a circuit as overdamped, critically damped, or underdamped from the discriminant $\Delta = R^2 - 4L/C$ .
Compute the undamped natural frequency $\omega_0 = 1/\sqrt{LC}$ , the damping ratio $\zeta = (R/2)\sqrt{C/L}$ , and the damped frequency $\omega_d = \omega_0\sqrt{1-\zeta^2}$ .
Track energy as it sloshes between the capacitor (electric field) and the inductor (magnetic field), with the resistor monotonically converting it into heat.
Implement both an analytic solver and a numerical state-space solver in Python, and reproduce them with PyTorch's matrix exponential.

Why an RLC Circuit Is the Same Story

"Electromagnetism and mechanics share a hidden mathematical spine." The RLC circuit is the cleanest place to feel this.

In Section 22.06 we studied a mass on a spring with damping. We are about to write down the equation for a circuit made of an inductor, a resistor, and a capacitor in series. Letter by letter the two systems are different — one has mass and forces, the other has charge and voltage. Coefficient by coefficient they are identical. The same characteristic polynomial appears, the same three damping regimes appear, the same energy-exchange picture appears.

That is not a coincidence. It is a consequence of the universal structure of any linear oscillator: inertia + restoring force + dissipation. Once you can solve one, you can solve the other — which is why a single section of calculus unlocks both your physics and your electrical-engineering courses.

The equation we are about to solve

For a series circuit with inductance $L$ , resistance $R$ , capacitance $C$ , and capacitor charge $q(t)$ :

L\,\ddot q + R\,\dot q + \tfrac{1}{C}\,q = 0.

The loop current is $i(t) = \dot q(t)$ . Two initial conditions are required: $q(0)$ (initial capacitor charge) and $i(0) = \dot q(0)$ (initial current).

The Three Elements

Before we derive the ODE, let's make sure each element's voltage law is crystal clear. Voltage $v$ across an element and current $i$ through it are related by:

Element	Voltage law	Plain-English meaning
Resistor (R)	v_R = R · i	Ohm's law. The resistor instantly converts current into heat in proportion to its current.
Inductor (L)	v_L = L · di/dt	An inductor resists CHANGES in current. The faster current changes, the bigger the voltage it develops to oppose that change (Lenz's law).
Capacitor (C)	v_C = q / C	A capacitor stores charge q. Its voltage is just charge divided by capacitance. Note: i = dq/dt, so charge builds up by integrating current.

Why the inductor is the "mass" of a circuit

The inductor formula $v_L = L\,\frac{di}{dt}$ looks exactly like Newton's second law $F = m\,\frac{dv}{dt}$ — force opposes changes in velocity. Inductors are mechanical inertia for electrons. The capacitor, meanwhile, plays the role of the spring: it stores energy by being "displaced" from zero charge.

Deriving the RLC Equation from Kirchhoff

Kirchhoff's voltage law (KVL) says: the sum of the voltage drops around any closed loop is zero. For our series RLC loop with no external source:

v_L + v_R + v_C = 0.

Substitute each element's law:

L\,\frac{di}{dt} \;+\; R\,i \;+\; \frac{q}{C} \;=\; 0.

Now the trick: everything in terms of one variable. Because current is the rate of charge flow, $i = dq/dt$ , and $di/dt = d^2 q / dt^2$ . Substituting:

\boxed{\;L\,\ddot q \;+\; R\,\dot q \;+\; \tfrac{1}{C}\,q \;=\; 0\;}

This is a homogeneous second-order linear ODE with constant coefficients. Compare with the abstract form $a y'' + b y' + c y = 0$ from Section 22.01:

Abstract symbol	RLC meaning
a	L (inductance)
b	R (resistance)
c	1/C (inverse capacitance, sometimes called 'elastance')
y(t)	q(t) (capacitor charge)
y'(t)	i(t) (loop current)

Characteristic Equation & Roots

Try $q(t) = e^{rt}$ . Then $\dot q = r e^{rt}$ and $\ddot q = r^2 e^{rt}$ . Substituting and dividing out $e^{rt}$ :

L\,r^2 + R\,r + \tfrac{1}{C} = 0,

a plain quadratic in $r$ with discriminant

\Delta \;=\; R^2 - \frac{4L}{C}.

The roots are

r_{1,2} = \frac{-R \pm \sqrt{R^2 - 4L/C}}{2L}.

Sign of $\Delta$ tells you everything. The boundary $\Delta = 0$ happens when

R \;=\; R_{\text{crit}} \;=\; 2\sqrt{L/C}.

Above $R_{\text{crit}}$ , the circuit is too lossy to oscillate — it just slumps back to zero charge. Below it, the circuit rings.

The two engineering parameters

Engineers almost never quote $(L, R, C)$ directly. They use:

$\omega_0 = 1/\sqrt{LC}$ — the undamped natural frequency. The frequency the circuit would oscillate at if R = 0.
$\zeta = \tfrac{R}{2}\sqrt{C/L} = \tfrac{R}{2L\omega_0}$ — the damping ratio. Dimensionless.

In these variables the ODE becomes $\ddot q + 2\zeta\omega_0\,\dot q + \omega_0^2\,q = 0$ , which is the universal form shared by every weakly-damped second-order system in nature.

The Three Damping Regimes

Regime	Condition	Roots	q(t)
Overdamped	ζ > 1 (R > 2√(L/C))	Two distinct negative reals r₁, r₂	C₁ e^(r₁t) + C₂ e^(r₂t) — no oscillation, slow return
Critically damped	ζ = 1 (R = 2√(L/C))	One repeated real root r = −R/(2L)	(C₁ + C₂ t) e^(rt) — fastest non-oscillating return
Underdamped	0 < ζ < 1 (R < 2√(L/C))	α ± iβ with α = −R/(2L), β = √(−Δ)/(2L)	e^(αt) (C₁ cos βt + C₂ sin βt) — damped sinusoid

In the underdamped case the angular oscillation frequency is $\omega_d = \beta = \omega_0\sqrt{1-\zeta^2}$ . This is lower than the undamped natural frequency $\omega_0$ — damping not only kills the amplitude, it slightly slows the ringing.

The s-plane picture

Engineers plot the roots in the complex s-plane. Left half-plane ⇒ stable (decaying). Right half-plane ⇒ unstable (growing). Imaginary axis ⇒ pure (undamped) oscillation. For our passive RLC with $R, L, C > 0$ , the roots are always in the left half-plane — passive circuits cannot go unstable on their own.

Mechanical–Electrical Analogy

Match the coefficients of $m\ddot x + c\dot x + kx = 0$ against $L\ddot q + R\dot q + (1/C)q = 0$ :

Mechanical	Symbol	Electrical analog	Symbol
Mass (inertia)	m	Inductance (electromagnetic inertia)	L
Damping coefficient	c	Resistance	R
Spring stiffness	k	Inverse capacitance (elastance)	1/C
Displacement	x	Capacitor charge	q
Velocity	ẋ	Current	i = q̇
Applied force	F(t)	Source voltage	v_s(t)
Kinetic energy	½ m ẋ²	Magnetic energy in inductor	½ L i²
Spring potential energy	½ k x²	Electric energy in capacitor	q²/(2C)

This is not a metaphor — it is a literal mathematical equivalence. Any theorem you prove about one system maps directly onto the other. Electrical engineers exploit this so thoroughly that simulators like SPICE solve the SAME ODE solver backend that mechanical-engineering FEM codes do.

The Energy View: An LC Tank Sloshes, R Burns

Multiply the ODE through by $\dot q = i$ :

L\,\ddot q\,\dot q + R\,\dot q^2 + \tfrac{1}{C}\,q\,\dot q = 0.

The first and last terms are perfect derivatives: $L\,\ddot q\,\dot q = \tfrac{d}{dt}\bigl(\tfrac{1}{2}L\dot q^2\bigr)$ and $\tfrac{1}{C}q\,\dot q = \tfrac{d}{dt}\bigl(\tfrac{q^2}{2C}\bigr)$ . So:

\frac{d}{dt}\underbrace{\Bigl[\tfrac{1}{2}L\,i^2 + \tfrac{q^2}{2C}\Bigr]}_{\text{total stored energy } U(t)} \;=\; -R\,i^2.

Read this carefully — it's one of the most beautiful equations in elementary physics:

If R = 0 the right-hand side is zero, so $U(t)$ is conserved. Energy oscillates perfectly between the capacitor (electric field, term $q^2/(2C)$ ) and the inductor (magnetic field, term $\tfrac{1}{2}L i^2$ ). It sloshes back and forth forever.
If R > 0 the right-hand side is strictly negative whenever $i \neq 0$ — energy drains monotonically out of the LC tank as heat in the resistor.
The instantaneous power dissipated by the resistor is $P_R(t) = R\,i(t)^2$ — Joule heating. Integrating this gives the total heat produced; it equals the total initial energy in the long-time limit.

LC tanks are oscillators

A capacitor and an inductor with no resistance form a perfect oscillator at frequency $\omega_0 = 1/\sqrt{LC}$ . This is the principle behind radio tuners, crystal oscillators, and the quartz-clock crystal in every watch — choose $L$ and $C$ to put $\omega_0$ at the station you want to receive.

Worked Example (Step-by-Step)

Solve $q'' + 0.4\,q' + 4\,q = 0$ with $q(0) = 1\,\text{C},\; i(0) = q'(0) = 0\,\text{A}$ . These coefficients correspond to $L = 1\,\text{H}$ , $R = 0.4\,\Omega$ , $C = 0.25\,\text{F}$ (so $1/C = 4$ ).

Click to expand the full hand calculation

Step 1 — Read off L, R, C.

Compare with $L\ddot q + R\dot q + (1/C)q = 0$ : $L = 1,\; R = 0.4,\; 1/C = 4$ so $C = 0.25$ .

Step 2 — Characteristic equation.

r^2 + 0.4\,r + 4 = 0.

Step 3 — Discriminant.

$\Delta = (0.4)^2 - 4(1)(4) = 0.16 - 16 = -15.84$ . Negative ⇒ underdamped (complex-conjugate roots, the circuit will ring).

Step 4 — Critical resistance check.

$R_{\text{crit}} = 2\sqrt{L/C} = 2\sqrt{1/0.25} = 2\cdot 2 = 4$ . Our $R = 0.4$ is far below $R_{\text{crit}}$ , so this is a lightly-damped circuit.

Step 5 — Engineering parameters.

$\omega_0 = 1/\sqrt{LC} = 1/\sqrt{0.25} = 2$ rad/s. $\zeta = R/(2L\omega_0) = 0.4/(2\cdot 1\cdot 2) = 0.1$ (dimensionless).

Step 6 — Complex roots.

$\alpha = -R/(2L) = -0.4/2 = -0.2$ (decay rate). $\beta = \sqrt{-\Delta}/(2L) = \sqrt{15.84}/2 \approx 1.9900$ (damped frequency).

Cross-check with $\omega_d = \omega_0\sqrt{1-\zeta^2} = 2\sqrt{1-0.01} = 2\cdot 0.99499 \approx 1.9900$ ✓.

Step 7 — General solution.

q(t) = e^{-0.2\,t}\bigl(A_1\cos(1.99\,t) + A_2\sin(1.99\,t)\bigr).

Step 8 — Apply $q(0) = 1$ .

At $t = 0$ : $q(0) = A_1 \cdot 1 + A_2 \cdot 0 = A_1$ ⇒ $A_1 = 1$ .

Step 9 — Differentiate and apply $i(0) = 0$ .

Using the product rule on $q(t) = e^{\alpha t}(A_1\cos\beta t + A_2\sin\beta t)$ :

\dot q(t) = e^{\alpha t}\bigl[(\alpha A_1 + \beta A_2)\cos\beta t + (\alpha A_2 - \beta A_1)\sin\beta t\bigr].

At $t = 0$ : $\dot q(0) = \alpha A_1 + \beta A_2 = 0$ ⇒ $A_2 = -\alpha A_1 / \beta = 0.2 / 1.99 \approx 0.1005$ .

Step 10 — Particular solution.

\boxed{\,q(t) \;=\; e^{-0.2\,t}\bigl(\cos(1.99\,t) + 0.1005\,\sin(1.99\,t)\bigr)\,}

Step 11 — Sanity check.

$q(0) = 1\cdot(1 + 0) = 1$ ✓. $i(0) = \alpha\cdot 1 + \beta\cdot 0.1005 = -0.2 + 1.99\cdot 0.1005 \approx -0.2 + 0.2 = 0$ ✓.

Step 12 — A long-term value.

At $t = 2$ : $e^{-0.4} \approx 0.6703$ , $\cos(3.98) \approx -0.6663$ , $\sin(3.98) \approx -0.7457$ . So

q(2) \approx 0.6703\cdot(-0.6663 + 0.1005\cdot(-0.7457)) \approx 0.6703\cdot(-0.7412) \approx -0.4970.

The capacitor has swung past zero and now holds about half a coulomb of opposite polarity — the LC tank is ringing. Verify this against the interactive plot below (L = 1, R = 0.4, C = 0.25, q(0) = 1, i(0) = 0).

Step 13 — Energy check.

Initial energy: $U(0) = q(0)^2/(2C) + \tfrac{1}{2}L i(0)^2 = 1/(2\cdot 0.25) + 0 = 2\,\text{J}$ . At $t = 2$ , the envelope has shrunk by $e^{-0.4} \approx 0.6703$ , so total energy is roughly $U(2) \approx 0.6703^2\cdot 2 = 0.898\,\text{J}$ . The missing $\approx 1.1\,\text{J}$ has been dissipated as heat in $R$ .

Interactive RLC Explorer

Drag the sliders for $L$ , $R$ , $C$ , $q(0)$ , and $i(0)$ and watch the capacitor charge polarity flip, current dots flow around the loop, and the time series react in real time. The readouts show the damping regime, the natural frequency $\omega_0$ , the damping ratio $\zeta$ , and the damped frequency $\omega_d$ .

Loading interactive RLC circuit…

Things to try

Set $R = 0$ for a perfect LC tank. The charge and current swap roles forever; the envelope is a flat horizontal line.
Crank $R$ up to $2\sqrt{L/C}$ (the critical value). The charge slumps to zero in the fastest possible non-oscillating way.
Push $R$ beyond that. Now the circuit is overdamped — too sluggish to swing past zero.
Shrink $C$ . The natural frequency $\omega_0 = 1/\sqrt{LC}$ increases, so the circuit rings faster.

Interactive Energy Visualizer

Same circuit, different lens. The three bars on top show the instantaneous electric energy in the capacitor (red), the magnetic energy in the inductor (purple), and the cumulative heat dumped into the resistor (orange). The plot underneath traces all three over time. With $R = 0$ the dashed total-energy curve is a flat line; with $R > 0$ it decays as $e^{-Rt/L}$ .

Loading energy visualization…

Why the total decays at e^(-Rt/L), not e^(-Rt/(2L))

The amplitude of $q(t)$ decays as $e^{\alpha t} = e^{-Rt/(2L)}$ . But energy is quadratic in $q$ and in $i$ , so the energy envelope decays at double the rate: $(e^{-Rt/(2L)})^2 = e^{-Rt/L}$ . This is why oscillators are often characterised by their quality factor $Q = \omega_0 L / R$ — the number of radians the system rings through before its energy drops by a factor of $e$ .

Computation in Python

We will write two solvers: one analytic (using the closed-form solutions we just derived for all three regimes), and one numerical (using classical RK4 on the equivalent 2-D first-order system). On any reasonable time grid they should agree to ~10⁻⁸.

Analytic solver — closed form for all three regimes

Direct evaluation of q(t) once the characteristic roots are known

🐍rlc_analytic.py

Explanation(15)

Code(72)

1Import NumPy

NumPy gives us elementwise exp, cos, sin, sqrt over the whole time grid. The entire solution is one vectorized expression — no Python loop.

4Function header — five arguments

Three circuit constants (L, R, C), two initial conditions (q(0) and i(0) = q′(0)), and a time array t. Why TWO initial conditions? Because the ODE is second-order. Physically: at t=0 we must know both how much charge is on the capacitor AND how fast it is flowing.

25Discriminant of the characteristic polynomial

The characteristic equation is L·r² + R·r + 1/C = 0. Its discriminant is Δ = R² − 4L/C. SIGN of Δ classifies the entire qualitative behaviour of the circuit. We do NOT yet need the roots themselves — just the sign.

27Branch 1 — overdamped (Δ > 0)

Two distinct REAL roots r₁, r₂. Both are negative (because R > 0 and L > 0, both coefficients of r have the same sign), so each mode decays as a plain exponential. No oscillation. Physically: the resistor dissipates energy faster than the LC tank can swap it.

29Two roots from the quadratic formula

r₁,₂ = (−R ± √Δ) / (2L). Each root is a damping rate (1/seconds). The MORE negative root corresponds to the FASTER-decaying mode.

32Initial-condition system

C₁ + C₂ = q(0) and r₁·C₁ + r₂·C₂ = q′(0) = i(0). Two linear equations in two unknowns. Solving gives C₂ first, then C₁ = q₀ − C₂. The denominator r₂ − r₁ ≠ 0 because we are in the distinct-roots branch.

34Vectorized evaluation

NumPy applies exp elementwise over the time grid t, so q is a 1-D array the same length as t. The current i = q′ is the derivative of the same combination of exponentials — multiply each term by its root.

37Branch 2 — underdamped (Δ < 0)

Complex-conjugate roots r = α ± iβ. Real part α = −R/(2L) is the envelope decay rate. Imaginary part β = √(−Δ)/(2L) is the damped angular frequency. We never write complex numbers in the final answer — Euler's formula turns the complex exponential into a REAL damped sinusoid.

39Real and imaginary parts of the root

α controls how fast the amplitude shrinks (envelope e^(αt)). β controls how fast the phase rotates inside that envelope (cos(βt), sin(βt)). They are SEPARATE physical quantities even though they emerge from the same characteristic root.

41Initial-condition system, underdamped case

At t=0, cos(0)=1 and sin(0)=0, so q(0) = A₁. Differentiating the product gives q′(0) = α·A₁ + β·A₂, hence A₂ = (i₀ − α·q₀) / β. The β in the denominator is nonzero because we are in the complex-roots branch (Δ < 0 ⇒ β > 0).

43Damped sinusoid

The factor env = e^(αt) is the SHRINKING envelope. Inside, (A₁ cos βt + A₂ sin βt) is a pure sinusoid that you can also write as M·cos(βt − φ) for some amplitude M and phase φ. The two forms are equivalent.

44Current i(t) for the underdamped case

Differentiate the product e^(αt)·(A₁ cos + A₂ sin). The product rule produces a combination of α·(...) and β·(derivative inside). Notice that current is NOT simply the derivative of the cosine part with α=0 — the envelope contributes too. This is why current and charge are NOT exactly π/2 out of phase when R > 0.

49Branch 3 — critically damped (Δ = 0)

Exactly ONE repeated root r = −R/(2L). One exponential e^(rt) is not enough to span the 2-D solution space, so we need a second linearly independent solution: t·e^(rt). The general solution is (C₁ + C₂t)·e^(rt). This is the boundary between oscillating and non-oscillating behaviour — the FASTEST return to zero charge without overshoot. Power-supply designers tune toward this regime on purpose.

56Run the worked example

L=1 H, R=0.4 Ω, C=0.25 F. Discriminant: (0.4)² − 4·1·(1/0.25) = 0.16 − 16 = −15.84. Negative → underdamped. α = −0.2, β = √15.84/2 ≈ 1.99. So q(t) = e^(−0.2t)·(cos(1.99t) + 0.1005·sin(1.99t)).

60Print key diagnostics

ω₀ = 1/√(LC) = 2 rad/s is the UNDAMPED natural frequency. ζ = (R/2)·√(C/L) = 0.1 is the dimensionless damping ratio. Because ζ < 1, the circuit rings before settling — and the printed value q(t=2) ≈ −0.498 confirms that at t=2 s the capacitor has swung past zero and is holding nearly half a coulomb of OPPOSITE polarity.

57 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def solve_rlc(L, R, C, q0, i0, t):
5    """
6    Solve the homogeneous series RLC equation
7        L q'' + R q' + (1/C) q = 0
8    for the capacitor charge q(t).
9
10    Parameters
11    ----------
12    L, R, C : float
13        Inductance (H), resistance (Ohm), capacitance (F).
14    q0 : float   q(0),  initial capacitor charge.
15    i0 : float   q'(0), initial loop current.
16    t  : np.ndarray, time grid.
17
18    Returns
19    -------
20    q : np.ndarray, charge on the capacitor at each t.
21    i : np.ndarray, current i(t) = q'(t).
22    info : dict, characteristic-equation diagnostics.
23    """
24    # characteristic polynomial:  L r^2 + R r + 1/C = 0
25    a, b, c = L, R, 1.0 / C
26    disc = b * b - 4 * a * c
27
28    if disc > 0:
29        # Overdamped: two distinct negative real roots
30        sD = np.sqrt(disc)
31        r1 = (-b + sD) / (2 * a)
32        r2 = (-b - sD) / (2 * a)
33        C2 = (i0 - r1 * q0) / (r2 - r1)
34        C1 = q0 - C2
35        q = C1 * np.exp(r1 * t) + C2 * np.exp(r2 * t)
36        i = C1 * r1 * np.exp(r1 * t) + C2 * r2 * np.exp(r2 * t)
37        return q, i, {"regime": "overdamped", "r1": r1, "r2": r2}
38
39    if disc < 0:
40        # Underdamped: damped sinusoid
41        alpha = -b / (2 * a)              # decay rate
42        beta  = np.sqrt(-disc) / (2 * a)  # damped angular frequency
43        A1 = q0
44        A2 = (i0 - alpha * q0) / beta
45        env = np.exp(alpha * t)
46        q = env * (A1 * np.cos(beta * t) + A2 * np.sin(beta * t))
47        i = env * (
48            (alpha * A1 + beta * A2) * np.cos(beta * t)
49            + (alpha * A2 - beta * A1) * np.sin(beta * t)
50        )
51        return q, i, {"regime": "underdamped", "alpha": alpha, "beta": beta}
52
53    # disc == 0: critically damped, single repeated root
54    r = -b / (2 * a)
55    C1 = q0
56    C2 = i0 - r * q0
57    q = (C1 + C2 * t) * np.exp(r * t)
58    i = (C2 + r * (C1 + C2 * t)) * np.exp(r * t)
59    return q, i, {"regime": "critical", "r": r}
60
61
62# ---- run the worked example ----
63L, R, C = 1.0, 0.4, 0.25
64q0, i0  = 1.0, 0.0
65t = np.linspace(0, 20, 800)
66q, i, info = solve_rlc(L, R, C, q0, i0, t)
67
68print(info)                              # {'regime': 'underdamped', 'alpha': -0.2, 'beta': 1.9899...}
69print("omega_0   =", 1 / np.sqrt(L * C)) # 2.0  (undamped natural frequency)
70print("zeta      =", (R / 2) * np.sqrt(C / L))  # 0.1
71print("q(t=2)    =", q[80])              # ~ -0.4983
72print("i(t=2)    =", i[80])

Numerical solver — RK4 on the 2-D state-space form

Any real-world RLC circuit you simulate (with diodes, op-amps, time- varying components, nonlinear loads) will require a numerical solver. The standard recipe converts the second-order ODE into a first-order vector ODE by stacking $Y = [q,\; i]^T$ :

\frac{d}{dt}\begin{bmatrix} q \\ i \end{bmatrix} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\tfrac{1}{LC} & -\tfrac{R}{L} \end{bmatrix}}_{A}\begin{bmatrix} q \\ i \end{bmatrix}.

The eigenvalues of this state matrix $A$ are exactly the characteristic roots $r_{1,2}$ . Engineers call this the state-space form of the circuit, and it is the universal input format for control-theory tools.

RK4 integration of the equivalent first-order vector ODE

🐍rlc_rk4.py

Explanation(8)

Code(45)

3RK4 step function

Classical Runge–Kutta 4 takes one explicit step of the ODE Y′ = f(t, Y). It evaluates the right-hand side at four carefully placed points (k1 at the start, k2 and k3 at the midpoint with slightly different y-values, k4 at the end) and combines them. Local error is O(h⁵), global error O(h⁴). This is the same routine inside scipy's solve_ivp default and MATLAB's ode45.

10Reduce 2nd-order ODE to a 2-D 1st-order system

Numerical solvers consume FIRST-order systems. Introduce i = q′. Then q′ = i and i′ = q″ = −(R/L)·i − (1/(LC))·q. Stacking [q, i] gives a 2-D linear ODE Y′ = A·Y with constant matrix A.

16The state matrix A

Top row [0, 1] encodes q′ = i. Bottom row encodes the original RLC equation, solved for i′. The EIGENVALUES of A are exactly the roots of the characteristic equation L·r² + R·r + 1/C = 0 — overdamped means two real eigenvalues, underdamped means a complex-conjugate pair.

21Right-hand side f(t, Y)

Because the ODE is autonomous (no explicit t-dependence) and linear, f(t, Y) is simply A·Y. The @ operator is NumPy matrix–vector multiplication.

23Initial state

The two initial conditions become the two components of Y₀ = [q₀, i₀]. Both are needed — even if i₀ = 0 you must pass it explicitly.

27Time-stepping loop

Walk along the time grid, applying RK4 to march Y forward. At each step we only need the previous state and the step size h = t[n] − t[n−1]. This is how every general-purpose ODE solver works under the hood.

32Return q and i separately

We slice the (T, 2) state array into the two columns. Column 0 is charge q(t); column 1 is current i(t). For oscilloscope-style plotting you typically want both.

38Numerical-vs-analytic sanity check

Because we already have a closed form from solve_rlc, this is a perfect regression test. On a fine grid RK4 should match the exact solution to ~1e-8. If your RK4 implementation drifts noticeably from the closed form, the bug is in the integrator — not in the physics.

37 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_rlc(L, R, C, q0, i0, t):
12    """
13    Integrate the series RLC equation by converting to a first-order system
14        Y = [q, i]     Y' = A @ Y
15    with
16            [  0       1   ]
17        A = [ -1/(LC) -R/L ]
18    """
19    A = np.array([
20        [0.0,        1.0],
21        [-1.0/(L*C), -R/L],
22    ])
23
24    def f(_t, Y):
25        return A @ Y
26
27    Y = np.array([q0, i0], dtype=float)
28    out = np.zeros((len(t), 2))
29    out[0] = Y
30    for n in range(1, len(t)):
31        h = t[n] - t[n-1]
32        Y = rk4_step(f, Y, t[n-1], h)
33        out[n] = Y
34    return out[:, 0], out[:, 1]   # q(t), i(t)
35
36# Numerical run of the worked example
37t = np.linspace(0, 20, 800)
38q_num, i_num = simulate_rlc(L=1.0, R=0.4, C=0.25, q0=1.0, i0=0.0, t=t)
39
40# Compare against the closed-form solution computed earlier
41from numpy.testing import assert_allclose
42# (assume q_exact, i_exact were computed by solve_rlc)
43# assert_allclose(q_num, q_exact, atol=1e-6)
44print("q_num[80] =", q_num[80])   # ~ -0.4983
45print("i_num[80] =", i_num[80])

PyTorch View: State-Space and Matrix Exponential

Because $Y' = A\,Y$ is linear with constant $A$ , an extraordinary fact applies: the exact solution is

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

where $e^{At}$ is the matrix exponential. PyTorch ships $\texttt{torch.linalg.matrix\_exp}$ , which evaluates this to machine precision. No numerical-integrator error.

Closed-form RLC solution via the matrix exponential

🐍rlc_matrix_exp.py

Explanation(7)

Code(34)

1Import PyTorch

torch.linalg.matrix_exp gives us the dense MATRIX exponential. For a LINEAR constant-coefficient system Y′ = A·Y, the EXACT solution is Y(t) = e^(A·t)·Y(0). No numerical integrator needed.

10Build the state matrix A

Same companion-matrix form as in the RK4 version. Switching to PyTorch lets us run on GPU, autodiff through L, R, C, q₀, i₀, and fit any of those parameters to experimental traces.

14Initial state Y(0)

Stack the two initial conditions into a single 2-vector. The whole RLC problem reduces to one linear-algebra question: 'starting from Y(0), where does the matrix flow e^(At) take us?'

17Matrix exponential at every time step

For each t_i we form the 2×2 matrix A·t_i, exponentiate it, and multiply by Y(0). The result is Y(t_i) = [q(t_i), i(t_i)]. The matrix exponential is to matrices what the scalar exp is to numbers — and it is the natural way to solve any constant-coefficient linear ODE.

21Slice columns into q and i

Y has shape (T, 2). Column 0 holds the charge trajectory; column 1 holds the current trajectory. Both are float64 for accuracy.

27Energy book-keeping

U_C = q²/(2C) is the electric energy stored in the capacitor. U_L = ½·L·i² is the magnetic energy in the inductor. With R > 0, U_C + U_L decays monotonically — every joule eventually leaves the circuit as heat in the resistor. With R = 0 it would be perfectly conserved (a lossless LC tank).

29Long-time check

At t = 20 s in our example (α = −0.2) the envelope has shrunk by e^(−4) ≈ 0.0183, so U has shrunk by e^(−8) ≈ 0.000335. The number printed is essentially 0 within float64 — the resistor has consumed all the initial energy.

27 lines without explanation

1import torch
2
3def rlc_state_trajectory(L, R, C, q0, i0, t):
4    """
5    Exact RLC solution via the matrix exponential.
6
7    Write the ODE as
8        d/dt [q, i]^T = A [q, i]^T,
9    then the exact closed form is
10        [q(t), i(t)]^T = exp(A * t) @ [q(0), i(0)]^T.
11    """
12    A = torch.tensor([
13        [0.0,                   1.0],
14        [-1.0 / (L * C),  -R / L],
15    ], dtype=torch.float64)
16    Y0 = torch.tensor([q0, i0], dtype=torch.float64)
17
18    # Vectorized matrix exponential, one per time step
19    Y = torch.stack([
20        torch.linalg.matrix_exp(A * ti) @ Y0
21        for ti in t
22    ])
23    return Y[:, 0], Y[:, 1]   # q(t), i(t)
24
25t = torch.linspace(0, 20, 800, dtype=torch.float64)
26q, i = rlc_state_trajectory(L=1.0, R=0.4, C=0.25, q0=1.0, i0=0.0, t=t)
27print("q(t=2)  =", q[80].item())   # ~ -0.4983
28print("i(t=2)  =", i[80].item())
29
30# Energy book-keeping
31U_C = q ** 2 / (2 * 0.25)
32U_L = 0.5 * 1.0 * i ** 2
33print("max |U_C + U_L| =", (U_C + U_L).max().item())
34print("U at t=20       =", (U_C[-1] + U_L[-1]).item())  # essentially 0 — all converted to heat

Why PyTorch and not just NumPy?

Three reasons. First, autodiff: with PyTorch you can differentiate $q(t)$ through $L, R, C$ and fit them to experimental oscilloscope traces by gradient descent. Second, GPU acceleration: batched RLC simulations (e.g. parameter sweeps for a circuit-design optimization) parallelize trivially. Third, this is the building block of linear neural ODEs — a learnable matrix $A$ whose eigenvalues encode the continuous-time dynamics of a network.

Where This Shows Up

📻 Radio & signal processing

Tuned LC tanks select a single AM/FM station by setting $\omega_0$ to the carrier frequency.
Band-pass and notch filters in analog audio.
Crystal oscillators in every microcontroller clock.

⚡ Power electronics

Snubber circuits clamp transients on switching power supplies.
LC filters smooth the output of buck/boost converters.
EMI suppression chokes on USB and HDMI cables.

🩺 Medical & sensors

MRI gradient coils — heavily damped LRC stages.
Inductive proximity sensors and metal detectors.
Wireless power-transfer pickup coils.

🤖 Control & ML

State-space form is the input to every control-theory tool (PID tuning, LQR, Kalman filters).
Linear neural ODEs and SSMs (Mamba, S4) reuse the exact $e^{At}$ machinery.
Battery and motor models in EV control software.

Common Pitfalls

Don't confuse i(0) with q(0)

A series RLC needs both initial values. A common student mistake is to specify only the capacitor's initial charge and forget about the initial current. Without $i(0)$ the solution is under-determined.

Sign convention on q and i

We chose $i = dq/dt$ with a specific loop-traversal direction. Reversing the assumed direction of $i$ flips the sign of the cross-term but gives the same q(t). Stay consistent within a single analysis.

The discriminant is R² − 4L/C, NOT R² − 4LC

Because the coefficient of $q$ is $1/C$ , the discriminant of the characteristic polynomial picks up the $1/C$ : $R^2 - 4L/C$ , equivalently $R^2 - 4L\cdot(1/C)$ . Writing $R^2 - 4LC$ is dimensionally wrong.

Parallel RLC ≠ series RLC

A parallel RLC (R, L, C all connected to the same two nodes) obeys a similar second-order ODE, but the roles of $R$ and $G = 1/R$ swap, and the natural variable is voltage $v(t)$ rather than charge. The shape of the math is identical; the coefficient mapping is different.

Summary

Every series RLC circuit (no source) obeys the same second-order linear ODE. The recipe to analyze it is universal:

Write Kirchhoff's voltage law and substitute $i = dq/dt$ to get $L\ddot q + R\dot q + (1/C)q = 0$ .
Form the characteristic equation $Lr^2 + Rr + 1/C = 0$ and compute the discriminant $\Delta = R^2 - 4L/C$ .
Classify (overdamped / critical / underdamped) and write the appropriate general solution.
Apply the two initial conditions $q(0), i(0)$ to solve for $C_1, C_2$ (or $A_1, A_2$ ).
Compute engineering parameters $\omega_0 = 1/\sqrt{LC}$ , $\zeta = (R/2)\sqrt{C/L}$ , $\omega_d = \omega_0\sqrt{1-\zeta^2}$ .
Track energy via $U = \tfrac{1}{2}Li^2 + q^2/(2C)$ ; $dU/dt = -R i^2$ .

Master equation at a glance

Quantity	Mechanical (Sec 22.06)	Electrical (this section)
ODE	m ẍ + c ẋ + k x = 0	L q̈ + R q̇ + (1/C) q = 0
Inertia coefficient	m	L
Damping coefficient	c	R
Restoring coefficient	k	1/C
Natural frequency	ω₀ = √(k/m)	ω₀ = 1/√(LC)
Damping ratio	ζ = c/(2√(mk))	ζ = (R/2)√(C/L)
Critical condition	c² = 4mk	R² = 4L/C (R_crit = 2√(L/C))
Energy	½ m ẋ² + ½ k x²	½ L i² + q²/(2C)
Dissipated power	c ẋ²	R i²

The slogan

"The inductor is the mass, the resistor is the friction, the capacitor is the spring."

Coming next: Chapter 22 wraps up here. With forced oscillations, resonance, mechanical vibrations, and RLC circuits behind us, you now own the entire theory of second-order linear ODEs with constant coefficients — the workhorse equation of classical physics, signal processing, and continuous-time control. Chapter 23 lifts the assumption that the coefficients are constant.