Simulating Electrical Systems

Electrical circuits are among the most natural applications of ordinary differential equation (ODE) simulation. The voltages and currents in resistor-capacitor (RC) and resistor-inductor-capacitor (RLC) networks evolve continuously in time according to well-understood physical laws — Kirchhoff's voltage and current laws, combined with the constitutive relationships of each component. Because these governing equations are ODEs with known analytical solutions in simple cases, electrical circuits provide an ideal testbed: you can build a numerical simulator, run it, and then compare its output directly against an exact formula to verify correctness. This topic develops that workflow from first principles, covering the physics of RC and RLC circuits, the translation of those physics into state-space form suitable for a numerical integrator, frequency-domain analysis of circuit responses, and rigorous validation against analytical solutions.

The RC Circuit

The simplest dynamic electrical circuit consists of a resistor R (in ohms, Ω) in series with a capacitor C (in farads, F), driven by a voltage source Vs. Applying Kirchhoff's voltage law (KVL) around the single loop — the source voltage must equal the sum of the voltage drops across the resistor and the capacitor — gives:

Vs = V_R + Vc
Vs = R·I + Vc

The current through the capacitor is related to the rate of change of the capacitor voltage by I = C · dVc/dt. Substituting this into the KVL equation and rearranging produces the governing first-order ODE:

R·C·dVc/dt + Vc = Vs

This single equation completely describes the evolution of Vc over time once an initial condition is specified. The product R·C has units of seconds and is called the time constant, denoted τ (tau):

τ = R · C

The time constant is the most important single parameter of an RC circuit because it sets the characteristic timescale of every transient. A large τ means the circuit responds slowly; a small τ means it responds quickly. Physically, τ is the time it takes the capacitor voltage to reach approximately 63.2% of its final value when charging from zero — a direct consequence of the exponential form of the solution.

For the specific case of a step input (Vs suddenly jumps from 0 to a constant value Vs at time t = 0) with the capacitor initially uncharged (Vc(0) = 0), the ODE has an exact analytical solution:

Vc(t) = Vs · (1 - e^(-t/τ))

This exponential charging curve is fundamental. At t = τ the voltage has reached Vs·(1 - e⁻¹) ≈ 0.632·Vs. At t = 2τ it has reached about 86.5% of Vs. After five time constants (t = 5τ) the capacitor voltage equals Vs·(1 - e⁻⁵) ≈ 0.9933·Vs, which is within 1% of the source voltage. This is why engineers use five time constants as the practical definition of "fully charged": the remaining error is less than 1%.

For a discharging scenario — Vs = 0 and initial condition Vc(0) = V₀ — the solution is instead:

Vc(t) = V₀ · e^(-t/τ)

The capacitor exponentially decays toward zero. Again, after five time constants the stored charge is less than 1% of its original value.

As a concrete numerical example, consider R = 1000 Ω (1 kΩ) and C = 100 × 10⁻⁶ F (100 µF). Then τ = 1000 × 100 × 10⁻⁶ = 0.1 s. With a step voltage of Vs = 5 V, the capacitor reaches about 3.16 V at t = 0.1 s and is considered fully charged by t = 0.5 s. These numbers are easy to verify with a voltmeter and stopwatch, making the RC circuit a perfect calibration case for any numerical simulator.

The RLC Circuit

Adding an inductor L (in henries, H) in series with the resistor and capacitor creates an RLC circuit. Inductors resist changes in current according to V_L = L · dI/dt. Applying KVL to the series RLC loop:

Vs = V_R + V_L + Vc
Vs = R·I + L·dI/dt + Vc

Because Vc relates to I through I = C · dVc/dt, differentiating the KVL equation with respect to time and substituting gives a second-order ODE in the loop current I:

L·d²I/dt² + R·dI/dt + I/C = dVs/dt

This equation has two roots that determine the character of the transient response. Defining two key parameters — the resonant frequency ω₀ and the damping factor α:

ω₀ = 1 / √(L·C)       [rad/s]
α  = R / (2·L)         [rad/s or s⁻¹]

The ratio of these quantities determines the circuit's damping regime. The characteristic equation of the homogeneous ODE is s² + 2α·s + ω₀² = 0, with roots s = −α ± √(α² − ω₀²). Three regimes arise:

• Overdamped (α > ω₀): Two distinct real negative roots. The current decays exponentially without oscillation — the circuit is sluggish.
• Critically damped (α = ω₀): A repeated real root. The current returns to steady state as fast as possible without overshooting — the fastest non-oscillatory response.
• Underdamped (α < ω₀): Two complex conjugate roots with negative real parts. The current oscillates at the damped natural frequency ω_d = √(ω₀² − α²) while its envelope decays exponentially. This is the oscillating transient behavior often seen in real circuits when R is small.

The underdamped case produces the characteristic ringing waveform. For example, with L = 10 mH, C = 1 µF, and R = 10 Ω: ω₀ = 1/√(0.01 × 10⁻⁶) = 10,000 rad/s (about 1592 Hz), and α = 10/(2 × 0.01) = 500 s⁻¹. Since α < ω₀, the circuit is underdamped and will ring at ω_d = √(10000² − 500²) ≈ 9987 rad/s.

The quality factor Q is a dimensionless number that quantifies how sharply the circuit resonates in the frequency domain:

Q = ω₀·L / R

A high Q (much greater than 1) means the circuit is very selective — it strongly amplifies signals near ω₀ and heavily attenuates all others. A low Q (less than 0.5) indicates an overdamped circuit with no resonance peak. At the borderline Q = 0.5 the circuit is critically damped. For the example above, Q = 10000 × 0.01 / 10 = 10, indicating a sharp resonance — exactly what you'd expect for a low-loss LC tank circuit.

Implementing the Electrical Simulation

To simulate either circuit numerically, the governing ODE(s) must be rewritten in state-space form: a vector of state variables whose time derivatives are expressed as explicit functions of the current state and time. A standard integrator such as scipy.integrate.solve_ivp (or a hand-coded RK4) then advances this state vector forward in time.

For the RC circuit, the state vector has a single element — the capacitor voltage Vc. The derivative function is obtained by rearranging the governing ODE:

state = [Vc]
dVc/dt = (Vs(t) - Vc) / (R * C)

In Python this becomes a function that accepts time t and state vector y, and returns the derivative:

def rc_deriv(t, y, R, C, Vs_func):
    Vc = y[0]
    Vs = Vs_func(t)
    dVc_dt = (Vs - Vc) / (R * C)
    return [dVc_dt]

For the RLC circuit, the second-order ODE must be split into two coupled first-order ODEs by introducing a second state variable — the derivative of current. Define the state as [I, dI/dt]:

state = [I, dI_dt]

d(I)/dt     = dI_dt
d(dI_dt)/dt = (dVs_dt - R*dI_dt - I/C) / L

The source term requires the time derivative of Vs. For a step input this is zero (after the initial jump); for a sinusoidal input Vs(t) = A·sin(ωt), the derivative is A·ω·cos(ωt). In code:

def rlc_deriv(t, y, R, L, C, dVs_func):
    I, dI_dt = y
    dVs = dVs_func(t)
    d2I_dt2 = (dVs - R * dI_dt - I / C) / L
    return [dI_dt, d2I_dt2]

Good simulation practice requires defining all physical component values as named constants at the top of the script, never embedding them as magic numbers inside functions. This makes the code readable, easy to modify, and less prone to errors from using inconsistent values in different places:

import numpy as np
from scipy.integrate import solve_ivp

# --- Component values ---
R = 1000.0       # Ohms
C = 100e-6       # Farads
L = 10e-3        # Henries (RLC only)
Vs_step = 5.0    # Volts (step amplitude)
tau = R * C      # Time constant

# --- Time array ---
t_start = 0.0
t_end   = 5 * tau
N       = 1000
t_eval  = np.linspace(t_start, t_end, N)

To fully characterize a circuit's behavior, two types of source inputs should be simulated. A step input (Vs suddenly switches from 0 to a constant value) reveals the transient response: how the circuit transitions from its initial state to a new steady state. A sinusoidal input Vs(t) = A·sin(2π·f·t) reveals the steady-state frequency response: how the circuit amplifies or attenuates signals of a given frequency after all transients have died away. By sweeping the frequency f across a range of values and recording the steady-state output amplitude, you can empirically construct the circuit's frequency response from the time-domain simulation.

Frequency Domain Analysis

Beyond time-domain transient plots, engineers analyze circuits in the frequency domain — decomposing the time-domain signal into its sinusoidal components and examining how the circuit modifies each one. Python's NumPy library provides the discrete Fourier transform (DFT) for this purpose.

np.fft.fft(signal) computes the DFT of an array signal, returning a complex array whose elements represent the amplitude and phase of each frequency component. The magnitude of each element gives the strength of that frequency in the signal:

import numpy as np

dt = t_eval[1] - t_eval[0]         # Sampling interval
Vc_fft = np.fft.fft(Vc_numerical)  # Complex spectrum
magnitude = np.abs(Vc_fft)         # Magnitude spectrum

np.fft.fftfreq(N, d=dt) returns the corresponding frequency axis in hertz, matching element-by-element with the FFT output array. Because the DFT of a real signal is symmetric, you typically keep only the first half of the spectrum (positive frequencies):

freqs = np.fft.fftfreq(N, d=dt)
# Keep only positive frequencies
pos_mask = freqs >= 0
freqs_pos = freqs[pos_mask]
magnitude_pos = magnitude[pos_mask]

The resulting magnitude spectrum reveals which frequency components are amplified or attenuated by the circuit. For an RC low-pass filter driven by a sinusoid near its cutoff frequency f_c = 1/(2πτ), the spectrum will show strong output below f_c and rapidly diminishing output above f_c. For an underdamped RLC circuit driven by a broadband signal (such as a short pulse), the spectrum will show a pronounced peak near the resonant frequency f₀ = ω₀/(2π).

The standard engineering representation of frequency response is the Bode plot, which displays two quantities as functions of frequency (on a logarithmic frequency axis):

• Magnitude (dB): 20 · log₁₀(|H(f)|), where H(f) is the ratio of output to input amplitude at frequency f. A value of 0 dB means no amplification or attenuation; −3 dB corresponds to the half-power (cutoff) point; −20 dB per decade is the slope of a first-order roll-off.
• Phase (degrees or radians): The phase shift between input and output sinusoids at each frequency.

You can construct an empirical Bode plot by running sinusoidal simulations at many frequencies, measuring the steady-state amplitude ratio and phase shift at each, and plotting the results. Alternatively, for linear circuits the theoretical transfer function can be evaluated analytically. For the RC low-pass filter, H(f) = 1 / (1 + j·2π·f·τ), giving |H| = 1/√(1 + (2πfτ)²) and phase = −arctan(2πfτ). At f = f_c = 1/(2πτ), |H| = 1/√2 ≈ 0.707, which is −3 dB, confirming the half-power cutoff.

import numpy as np
import matplotlib.pyplot as plt

f_range = np.logspace(0, 5, 500)     # 1 Hz to 100 kHz
H_mag = 1.0 / np.sqrt(1 + (2*np.pi*f_range*tau)**2)
H_phase = -np.arctan(2*np.pi*f_range*tau) * 180/np.pi
H_dB = 20 * np.log10(H_mag)

fig, (ax1, ax2) = plt.subplots(2, 1)
ax1.semilogx(f_range, H_dB)
ax1.set_ylabel('Magnitude (dB)')
ax2.semilogx(f_range, H_phase)
ax2.set_ylabel('Phase (deg)')
ax2.set_xlabel('Frequency (Hz)')
plt.tight_layout()
plt.show()

Validation Against Analytical Solutions

One of the greatest advantages of simulating circuits with known analytical solutions is the ability to rigorously validate the numerical method. Validation should always be quantitative, not merely visual — a plot where two curves appear to overlap is not sufficient evidence that the error is small enough for your purpose.

The first step is to compute the analytical solution over exactly the same time array used for the numerical simulation. Using NumPy array operations, this is efficient and vectorized:

# RC circuit analytical solution (step input, zero initial condition)
Vc_analytical = Vs_step * (1 - np.exp(-t_eval / tau))

# Numerical solution from solve_ivp
sol = solve_ivp(rc_deriv, [t_start, t_end], [0.0],
                t_eval=t_eval, args=(R, C, lambda t: Vs_step),
                method='RK45', rtol=1e-8, atol=1e-10)
Vc_numerical = sol.y[0]

Plot both solutions on the same axes using distinct line styles so they are distinguishable even in black-and-white print. A common convention is to plot the analytical solution as a solid line and the numerical solution as a dashed line or with markers:

plt.figure()
plt.plot(t_eval, Vc_analytical, 'k-',  linewidth=2, label='Analytical')
plt.plot(t_eval, Vc_numerical,  'r--', linewidth=1.5, label='Numerical (RK45)')
plt.xlabel('Time (s)')
plt.ylabel('Capacitor Voltage (V)')
plt.legend()
plt.title('RC Circuit: Analytical vs Numerical')
plt.show()

Next, compute the absolute pointwise error at every time step:

error = np.abs(Vc_numerical - Vc_analytical)

plt.figure()
plt.semilogy(t_eval, error)
plt.xlabel('Time (s)')
plt.ylabel('Absolute Error (V)')
plt.title('Numerical Error vs Time')
plt.show()

Plotting the error on a logarithmic y-axis is often more informative than a linear scale because the error typically spans many orders of magnitude. Examining this plot tells you where in time the integrator struggles most. For the RC circuit with a step input, the error is usually largest near t = 0 where the solution has its steepest gradient, and decays toward machine precision as the solution flattens out.

Finally, report scalar validation metrics that give a single-number characterization of overall accuracy:

• Maximum absolute error (MAE): np.max(error) — the worst-case pointwise deviation. This is the relevant metric when you need a guaranteed bound on the error at any instant.
• Root-mean-square error (RMSE): np.sqrt(np.mean(error**2)) — the average error across the entire simulation. This is more representative of typical accuracy than the maximum, which may be due to a single outlier near a discontinuity.

max_error  = np.max(error)
rmse_error = np.sqrt(np.mean(error**2))
print(f"Maximum absolute error: {max_error:.6e} V")
print(f"Root-mean-square error: {rmse_error:.6e} V")

For a well-implemented RK45 integrator with tight tolerances (rtol=1e-8, atol=1e-10) applied to a smooth RC circuit problem, you should expect a maximum error on the order of 10⁻⁸ V or smaller — comparable to the specified tolerance. If the error is much larger, it indicates either a coding error in the derivative function, inconsistent component values between the numerical and analytical setups, or tolerances that are too loose. This kind of systematic validation, performed before trusting any simulation result, is standard engineering practice and ensures that the numerical tool is working correctly before it is applied to problems where no analytical solution exists.