Scientific Machine Learning:
Examples of Numerical Analysis
This notebook puts the methods of the previous notebook to work on four engineering problems: the operating point of a diode, the energy of a measured signal, the motion of a damped oscillator, and the propagation of a wave along a string. Each problem is solved from scratch in NumPy, so that the mechanics of root finding, integration and time stepping are seen on a concrete case before they are hidden inside a library call or inside a physics-informed loss.
Learning Objectives¶
After studying this notebook, you should be able to
Nonlinear root finding
- formulate an engineering operating-point problem as a nonlinear equation $f(x)=0$;
- apply the bisection method from a sign-changing interval, and explain why it cannot fail;
- apply the Newton-Raphson method using an analytical derivative, and recognize its rapid convergence once the estimate is close enough;
- apply the secant method without evaluating a derivative at all;
- compare the convergence behaviour of the three methods, and interpret residuals and stopping criteria.
Numerical integration
- express an engineering quantity as a definite integral and construct a uniform grid for it;
- apply the trapezoidal rule and Simpson's rule directly from their formulas;
- compare a numerical result against an analytical benchmark, and measure how the number of sample points affects the error;
- recognize that a higher-order rule is not automatically the more accurate one, and that the structure of the integrand decides.
ODE numerical solution
- convert a second-order equation of motion into first-order state-space form;
- implement the explicit Euler step and march it forward in time;
- write the exact solution of a damped forced oscillator and use it as a reference;
- predict the stability limit from the eigenvalues of the system matrix, and explain why the Euler solution converges below it and diverges above it.
1D wave equation
- discretize the wave equation with central differences in time and space, and explain why the scheme needs three time levels and how the first step is started;
- define the Courant number and derive the CFL condition $C\le1$ from the amplification factor;
- state the domain-of-dependence argument that gives the same condition physically;
- explain why the accuracy of this scheme improves as $C\to1$ and why it is exact at $C=1$;
- verify the second-order convergence of the scheme by refining the mesh.
Table of Contents¶
Nonlinear Root Finding Diode operating point
Numerical Integration Mean-square value of a periodic signal
ODE Numerical Solution Harmonic forced vibration of a 1-DOF system
1D Wave Equation Transverse vibration of a stretched string
Problem: Diode operating point
Many engineering models lead to a nonlinear algebraic equation
$$ f(x)=0. $$
In this lab, we solve a nonlinear equation that determines the operating point of a diode circuit. The problem is simple enough to visualize, but it also illustrates the essential ideas behind three fundamental root-finding methods:
- Bisection method
- Newton-Raphson method
- Secant method
The emphasis in this lab is on implementing the numerical algorithms directly and observing how the sequence of approximations converges to the root.
Consider a simple diode circuit with a voltage source $V_S$, a resistor $R$, and a diode.
The diode current is modeled by the Shockley equation
$$ I_D = I_S\left(e^{V_D/V_T}-1\right), $$
where $V_D$ is the diode voltage, $I_S$ is the saturation current, and $V_T$ is the thermal voltage.
The current through the resistor is
$$ I_R = \frac{V_S-V_D}{R}. $$
At the operating point, the diode current and resistor current are equal:
$$ I_D=I_R. $$
Therefore, the operating-point voltage is obtained by solving
$$ f(V_D) = I_S\left(e^{V_D/V_T}-1\right) - \frac{V_S-V_D}{R} =0. $$
The numerical task is to determine the value of $V_D$ that makes this equation zero.
Why is this a root-finding problem?
The diode equation contains the exponential term $e^{V_D/V_T}$. Therefore, the operating-point equation cannot be rearranged into a simple linear expression for $V_D$.
Instead, we define
$$ f(V_D) = I_S\left(e^{V_D/V_T}-1\right) - \frac{V_S-V_D}{R}, $$
and search for the voltage at which
$$ f(V_D)=0. $$
Each numerical method in this lab generates a sequence of voltage estimates
$$ V_D^{(0)},V_D^{(1)},V_D^{(2)},\ldots $$
that should approach the same physical operating point.
The circuit and numerical parameters are taken directly from the lab problem.
| Quantity | Symbol | Value |
|---|---|---|
| Supply voltage | $V_S$ | $5.0$ V |
| Resistance | $R$ | $100$ $\Omega$ |
| Saturation current | $I_S$ | $10^{-12}$ A |
| Thermal voltage | $V_T$ | $0.026$ V |
| Initial interval | $[V_a,V_b]$ | $[0.60,0.67]$ V |
| Current residual tolerance | $\varepsilon_f$ | $10^{-8}$ A |
| Voltage tolerance | $\varepsilon_V$ | $10^{-8}$ V |
| Maximum updates | 50 |
For Newton-Raphson, the upper end of the initial interval, $V_D=0.67$ V, is used as the initial estimate.
For the Secant method, the two endpoints $0.60$ V and $0.67$ V are used as the initial estimates.
import numpy as np
import matplotlib.pyplot as plt
The four numbers below define the physical problem. Changing any of them moves the operating point, but none of the algorithms has to be modified.
# Circuit parameters
VS = 5.0 # supply voltage [V]
R = 100.0 # resistance [ohm]
IS = 1e-12 # saturation current [A]
VT = 0.026 # thermal voltage [V]
These control the algorithms rather than the physics: where the search starts, how small the residual and the voltage change must become before the iteration is accepted, and how many updates are allowed before the loop gives up.
# Root-finding parameters
Va = 0.60 # left end of the starting interval [V]
Vb = 0.67 # right end of the starting interval [V]
tol_f = 1e-8 # residual tolerance [A]
tol_V = 1e-8 # voltage tolerance [V]
nmax = 50 # maximum number of updates
All three methods evaluate the same function repeatedly, so it is defined once here and reused. The residual is the difference between the diode current and the resistor current,
$$ f(V_D) = \underbrace{I_S\left(e^{V_D/V_T}-1\right)}_{\text{diode current}} - \underbrace{\frac{V_S-V_D}{R}}_{\text{resistor current}} , $$
and it has the units of current. A value of $f$ close to zero means that the two currents nearly balance, which is what the operating point requires.
def f(V):
# current mismatch f(V) in amperes; zero at the operating point
ID = IS*(np.exp(V/VT) - 1)
IR = (VS - V)/R
return ID - IR
Derivative. The Newton-Raphson method needs the slope of the residual. Differentiating the expression above with respect to $V_D$ gives
$$ f'(V_D) = \frac{I_S}{V_T}e^{V_D/V_T} + \frac{1}{R}. $$
Both terms are positive, so $f$ is strictly increasing and the root is unique. This is a useful property: it guarantees that the sign-change test below identifies the operating point without ambiguity.
def df(V):
# derivative f'(V) in amperes per volt
dID = (IS/VT)*np.exp(V/VT)
dIR = 1/R # the resistor branch enters with a plus sign
return dID + dIR
The Bisection method requires an interval in which the function changes sign.
For the initial interval $[V_a,V_b]$, we therefore check
$$ f(V_a)f(V_b)<0. $$
This condition indicates that the interval brackets a root, provided the function is continuous over the interval.
fa = f(Va)
fb = f(Vb)
print('f(%.2f) = %.6e A' % (Va, fa))
print('f(%.2f) = %.6e A' % (Vb, fb))
print('Sign change:', fa*fb < 0)
The Bisection method maintains an interval $[a_n,b_n]$ that brackets the root.
At each iteration, compute the midpoint
$$ c_n=\frac{a_n+b_n}{2}. $$
Then evaluate $f(c_n)$.
If
$$ f(a_n)f(c_n)<0, $$
the root remains in $[a_n,c_n]$.
Otherwise, the root remains in $[c_n,b_n]$.
The interval is therefore reduced by approximately one half at every update.
In this lab, convergence requires both
$$ |f(c_n)|\le \varepsilon_f $$
and
$$ \frac{b_n-a_n}{2}\le \varepsilon_V. $$
The first condition measures the equation residual, while the second measures the uncertainty in the root location.
The residual at the left endpoint is stored once, so that the sign test inside the loop does not recompute it, and the midpoints and residuals are recorded as the iteration proceeds for the convergence plots of Sections 1.9 and 1.10.
# Bisection: initialization
a = Va
b = Vb
fa = f(a)
V_bis = [] # history of midpoints
f_bis = [] # history of |f| at those midpoints
# Bisection: iteration
for update in range(nmax):
Vm = (a + b)/2
fm = f(Vm)
V_bis.append(Vm)
f_bis.append(abs(fm))
# stop when the interval is narrow enough and the residual is small enough
interval_ok = (b - a)/2 <= tol_V
residual_ok = abs(fm) <= tol_f
if fm == 0 or (interval_ok and residual_ok):
break
# keep the half that still contains a sign change
if fa*fm < 0:
b = Vm
else:
a = Vm
fa = fm
V_root_bis = Vm
# Bisection: result
print('Root = %.10f V' % V_root_bis)
print('|f(root)| = %.3e A' % f_bis[-1])
print('Number of updates = %d' % len(V_bis))
Bisection result
The Bisection method gives
$$ V_D \approx 0.6369744360\ \text{V}. $$
The final residual is approximately
$$ |f(V_D)|\approx 3.82\times10^{-9}\ \text{A}. $$
The method requires 23 updates for the specified stopping criteria.
The important feature of Bisection is reliability: once a valid sign-changing interval has been found, the interval is systematically reduced and the root remains bracketed.
The Newton-Raphson method uses the tangent line to the nonlinear function.
Starting from $V_n$, the next estimate is
$$ V_{n+1} = V_n-\frac{f(V_n)}{f'(V_n)}. $$
For the diode operating-point equation,
$$ f(V_D) = I_S\left(e^{V_D/V_T}-1\right) - \frac{V_S-V_D}{R}, $$
the derivative is
$$ f'(V_D) = \frac{I_S}{V_T}e^{V_D/V_T} + \frac{1}{R}. $$
In this lab, the initial estimate is
$$ V_D^{(0)}=0.67\ \text{V}. $$
The method stops when both the residual and the change in voltage are sufficiently small.
Newton-Raphson needs a single starting voltage, taken here as the upper end of the interval used for bisection. The update count below excludes that initial estimate, so that it can be compared directly with the number of bisection updates.
# Newton-Raphson: initialization
V = Vb
V_nr = [V] # history of estimates
f_nr = [] # history of |f| at those estimates
# Newton-Raphson: iteration
for update in range(nmax):
V_next = V - f(V)/df(V) # tangent line crosses zero here
f_next = f(V_next)
V_nr.append(V_next)
f_nr.append(abs(f_next))
# stop when the estimate has settled and the residual is small enough
step_ok = abs(V_next - V) <= tol_V
residual_ok = abs(f_next) <= tol_f
if step_ok and residual_ok:
break
V = V_next
V_root_nr = V_next
# Newton-Raphson: result
print('Root = %.10f V' % V_root_nr)
print('|f(root)| = %.3e A' % f_nr[-1])
print('Number of updates = %d' % (len(V_nr) - 1))
Newton-Raphson result
The Newton-Raphson method gives
$$ V_D \approx 0.6369744338\ \text{V}. $$
The final residual is approximately
$$ |f(V_D)|\approx 5.55\times10^{-17}\ \text{A}. $$
Only 6 updates are required after the initial estimate.
This much faster convergence is the main advantage of Newton-Raphson in this example. However, unlike Bisection, Newton-Raphson does not maintain a root-bracketing interval at every step. Its performance therefore depends more strongly on the starting point and the local behavior of the function.
The Secant method avoids the explicit evaluation of $f'(V)$.
Instead, the derivative is approximated using two consecutive points:
$$ f'(V_n) \approx \frac{f(V_n)-f(V_{n-1})}{V_n-V_{n-1}}. $$
Substituting this approximation into the Newton update gives
$$ V_{n+1} = V_n - f(V_n) \frac{V_n-V_{n-1}} {f(V_n)-f(V_{n-1})}. $$
The Secant method therefore requires two initial estimates. In this lab,
$$ V_0=0.60\ \text{V}, \qquad V_1=0.67\ \text{V}. $$
Unlike Bisection, the Secant method does not explicitly preserve a sign-changing interval.
The secant method needs two starting voltages, because its slope is built from the two most recent points. Both are recorded in the history so that the plotted sequence begins where the algorithm begins, and both are excluded from the update count.
# Secant: initialization
V0 = Va
V1 = Vb
f0 = f(V0)
f1 = f(V1)
V_sec = [V0, V1] # history of estimates
f_sec = [abs(f0), abs(f1)] # history of |f| at those estimates
# Secant: iteration
for update in range(nmax):
# the line through the last two points crosses zero here
V_next = V1 - f1*(V1 - V0)/(f1 - f0)
f_next = f(V_next)
V_sec.append(V_next)
f_sec.append(abs(f_next))
step_ok = abs(V_next - V1) <= tol_V
residual_ok = abs(f_next) <= tol_f
if step_ok and residual_ok:
break
# shift the two points forward
V0 = V1
f0 = f1
V1 = V_next
f1 = f_next
V_root_sec = V_next
# Secant: result
print('Root = %.10f V' % V_root_sec)
print('|f(root)| = %.3e A' % f_sec[-1])
print('Number of updates = %d' % (len(V_sec) - 2))
Secant result
The Secant method gives
$$ V_D \approx 0.6369744338\ \text{V}. $$
The final residual is approximately
$$ |f(V_D)|\approx 7.05\times10^{-13}\ \text{A}. $$
The method requires 8 updates after the two initial estimates.
The Secant method reaches the solution much faster than Bisection in this example, while avoiding the need to derive and evaluate an analytical derivative.
All three methods converge to essentially the same physical operating point.
The small differences in the reported digits are caused by the different numerical sequences and stopping points. The important quantities for comparison are the root estimate, residual, and number of updates.
# Compare the three methods
n_bis = len(V_bis)
n_nr = len(V_nr) - 1
n_sec = len(V_sec) - 2
print('Method Root (V) |f(root)| (A) Updates')
print('-'*70)
print('Bisection %.10f %.3e %d' % (V_root_bis, f_bis[-1], n_bis))
print('Newton-Raphson %.10f %.3e %d' % (V_root_nr, f_nr[-1], n_nr))
print('Secant %.10f %.3e %d' % (V_root_sec, f_sec[-1], n_sec))
1.8.1. Summary¶
| Method | Root estimate (V) | Final residual (A) | Updates |
|---|---|---|---|
| Bisection | 0.6369744360 | $3.82\times10^{-9}$ | 23 |
| Newton-Raphson | 0.6369744338 | $5.55\times10^{-17}$ | 6 |
| Secant | 0.6369744338 | $7.05\times10^{-13}$ | 8 |
The three methods agree on the operating point:
$$ V_D \approx 0.636974434\ \text{V} $$
The comparison also illustrates a fundamental trade-off:
- Bisection: slower, but robust and bracket-preserving.
- Newton-Raphson: very fast near the root, but requires a derivative and a suitable initial estimate.
- Secant: fast and derivative-free, but does not provide the same bracketing guarantee as Bisection.
A useful way to understand a numerical method is to examine the sequence of estimates generated by the algorithm.
The following plot shows the estimated diode voltage versus update number. All three methods approach the same root, but they do so at different rates.
# Plot the voltage estimates against the update number
k_bis = range(1, len(V_bis) + 1)
k_nr = range(0, len(V_nr))
k_sec = range(0, len(V_sec))
plt.figure(figsize = (6, 4))
plt.plot(k_bis, V_bis, 'o-', label = 'Bisection')
plt.plot(k_nr, V_nr, 's-', label = 'Newton-Raphson')
plt.plot(k_sec, V_sec, '^-', label = 'Secant')
plt.axhline(0.6369744338, color = 'k', linestyle = '--', label = 'Reference root')
plt.xlabel('Iteration / update')
plt.ylabel('Diode voltage $V_D$ (V)')
plt.title('Convergence of the voltage estimates')
plt.grid(True)
plt.legend()
plt.show()
The residual
$$ |f(V_D)| $$
measures how closely a numerical estimate satisfies the original nonlinear equation.
A small residual is therefore an important indicator of convergence.
Because the residual can decrease over several orders of magnitude, a logarithmic vertical axis is useful for comparing the methods.
# Plot the residual against the update number, on a logarithmic axis
k_bis = range(1, len(f_bis) + 1)
k_nr = range(1, len(f_nr) + 1)
k_sec = range(1, len(f_sec) + 1)
plt.figure(figsize = (6, 4))
plt.semilogy(k_bis, f_bis, 'o-', label = 'Bisection')
plt.semilogy(k_nr, f_nr, 's-', label = 'Newton-Raphson')
plt.semilogy(k_sec, f_sec, '^-', label = 'Secant')
plt.axhline(tol_f, color = 'k', linestyle = '--', label = 'Residual tolerance')
plt.xlabel('Update')
plt.ylabel('$|f(V_D)|$ (A)')
plt.title('Residual convergence')
plt.grid(True)
plt.legend()
plt.show()
This lab demonstrates that different numerical algorithms can solve the same nonlinear engineering problem while exhibiting very different convergence behavior.
The Bisection method reduces the root-containing interval in a predictable way. Its main strength is robustness.
Newton-Raphson uses local derivative information. Once the estimate is sufficiently close to the root, the convergence can be very rapid.
The Secant method obtains similar local information without explicitly calculating the derivative. It uses two previous points to construct an approximate slope.
For this diode problem, all three methods converge to approximately
$$ V_D \approx 0.636974434\ \text{V}. $$
The numerical method should therefore be selected not only by the number of iterations, but also by the information available about the problem, the reliability required, and the quality of the initial estimates.
Problem: Mean-square value of a periodic signal
Numerical integration is widely used in engineering when an integral cannot be evaluated easily by hand, or when only sampled data are available.
In this lab, we compute the mean-square value of a periodic signal using two classical numerical integration methods:
- Trapezoidal Rule
- Simpson's Rule
The analytical result is used as a benchmark so that the numerical accuracy of the two methods can be evaluated directly.
Consider the periodic signal
$$ x(t)=6\cos(0.5t+1.5)\cos(2t-6), \qquad -2\pi \leq t \leq 2\pi $$
We want to calculate its mean-square value over one period:
$$ M_x = \frac{1}{T_0} \int_{-2\pi}^{2\pi}x(t)^2\,dt. $$
The integration interval has length
$$ T_0=2\pi-(-2\pi)=4\pi. $$
Therefore,
$$ M_x = \frac{1}{4\pi} \int_{-2\pi}^{2\pi}x(t)^2\,dt. $$
The corresponding RMS value is
$$ x_{\mathrm{RMS}}=\sqrt{M_x}. $$
For this signal, the analytical mean-square value is
$$ M_x=9, $$
so the analytical RMS value is
$$ x_{\mathrm{RMS}}=3. $$
These analytical values provide a reference against which the numerical results can be compared.
Where the analytical value comes from
The product of two cosines can be rewritten as a sum,
$$ x(t)=3\cos(2.5t-4.5)+3\cos(1.5t-7.5), $$
so the signal is the sum of two sinusoids, each of amplitude $3$. Over the interval $[-2\pi,2\pi]$ the first completes $5$ full cycles and the second completes $3$, and the cross term averages to zero because the two frequencies differ. Each sinusoid therefore contributes $3^2/2=4.5$ to the mean square, and
$$ M_x=4.5+4.5=9. $$
This decomposition will be useful again in Section 2.8, where the accuracy of the two rules turns out to depend on exactly which frequencies the signal contains.
Engineering interpretation
The mean-square value is an important quantity in signal processing and engineering because it measures the average squared magnitude of a signal.
For example, the RMS value of a voltage or current is directly related to average power in resistive systems.
In this lab, the main numerical task is therefore
$$ \int_{-2\pi}^{2\pi}x(t)^2\,dt. $$
Once this integral is approximated numerically, the mean-square and RMS values follow immediately.
We use uniformly spaced sample points over
$$ [-2\pi,\,2\pi]. $$
For Simpson's Rule, the number of sample points must be odd because the method operates on pairs of subintervals.
For the main calculation, we use
$$ N=31 $$
sample points.
The analytical benchmark is
$$ M_x=9. $$
# Integration interval
ta = -2*np.pi
tb = 2*np.pi
T0 = tb - ta # length of one period [s]
# Analytical benchmark
Mx_exact = 9.0
rms_exact = 3.0
def signal(t):
return 6*np.cos(0.5*t + 1.5)*np.cos(2*t - 6)
def integrand(t):
return signal(t)**2
The sample points are equally spaced and include both endpoints, so $N$ points give $N-1$ subintervals of width $h=(t_{\text{end}}-t_{\text{start}})/(N-1)$. With $N=31$ there are $30$ subintervals, an even number, so Simpson's Rule applies directly.
# Simpson's Rule requires an odd number of sample points
N = 31
# Uniform time grid
t = np.linspace(ta, tb, N)
x = signal(t)
x2 = integrand(t)
h = t[1] - t[0]
print('Number of sample points = %d' % N)
print('Number of subintervals = %d' % (N - 1))
print('Grid spacing h = %.6f' % h)
print('Integration interval = [%.4f, %.4f]' % (ta, tb))
print('Period T0 = %.6f' % T0)
Before performing numerical integration, it is useful to inspect the signal and the quantity that is actually being integrated.
The original signal is
$$ x(t), $$
while the integral uses
$$ x(t)^2. $$
Because squaring removes the sign of the signal, the integrand is always nonnegative.
t_dense = np.linspace(ta, tb, 500)
plt.figure(figsize = (6, 6))
plt.subplot(2, 1, 1)
plt.plot(t_dense, signal(t_dense), 'k', linewidth = 1, label = 'exact')
plt.plot(t, x, 'o', markersize = 4, label = '%d samples' % N)
plt.ylabel('x(t)')
plt.title('Signal x(t)')
plt.grid(True)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(t_dense, integrand(t_dense), 'k', linewidth = 1, label = 'exact')
plt.plot(t, x2, 'o', markersize = 4, label = '%d samples' % N)
plt.xlabel('Time t')
plt.ylabel('$x(t)^2$')
plt.title('Integrand $x(t)^2$')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
The integrand oscillates considerably faster than the signal, because squaring doubles the frequencies present. Thirty subintervals resolve it, but not generously, which is worth keeping in mind when the results below turn out to be far more accurate than that appearance suggests.
The Trapezoidal Rule approximates the area under a curve by replacing each small interval with a trapezoid.
For uniformly spaced points,
$$ t_i=t_0+ih, $$
the integral is approximated by
$$ I_T = h \left[ \frac{f_0}{2} + \sum_{i=1}^{N-2}f_i + \frac{f_{N-1}}{2} \right]. $$
For this problem,
$$ f_i=x(t_i)^2. $$
The mean-square value is then
$$ M_x^{(T)} = \frac{I_T}{T_0}. $$
The key idea is simple: interior points receive full weight, while the two endpoints receive half weight.
def trapz_Mx(N):
# mean-square value computed with the trapezoidal rule on N points
t = np.linspace(ta, tb, N)
y = integrand(t)
h = t[1] - t[0]
integral = h*(y[0]/2 + np.sum(y[1:-1]) + y[-1]/2)
return integral/T0
Mx_trap = trapz_Mx(N)
rms_trap = np.sqrt(Mx_trap)
err_trap = abs(Mx_trap - Mx_exact)
print('Trapezoidal mean-square value = %.10f' % Mx_trap)
print('Trapezoidal RMS value = %.10f' % rms_trap)
print('Absolute error = %.6e' % err_trap)
Trapezoidal result
For $N=31$ uniformly spaced points, the Trapezoidal Rule reproduces the analytical value to the last stored digit,
$$ M_x^{(T)}=9.000000000000, \qquad |M_x^{(T)}-9|\approx5\times10^{-15}. $$
That is round-off, not discretization error: the rule has effectively made no approximation at all. This is a surprising outcome for a method whose error is normally $O(h^2)$, and with $h\approx0.42$ on a rapidly oscillating integrand a visible error would have been expected. Section 2.8 explains why it does not appear.
Simpson's Rule improves the approximation by fitting a quadratic polynomial over pairs of adjacent subintervals.
For uniformly spaced points, the composite Simpson formula is
$$ I_S = \frac{h}{3} \left[ f_0 + 4\sum f_{\mathrm{odd}} + 2\sum f_{\mathrm{even}} + f_{N-1} \right]. $$
The weighting pattern is therefore
$$ 1,\ 4,\ 2,\ 4,\ 2,\ \ldots,\ 4,\ 1. $$
For Simpson's Rule, the number of subintervals must be even. Consequently, the number of sample points must be odd:
$$ N=2m+1. $$
With $N=31$, there are 30 subintervals, so Simpson's Rule can be applied directly.
The loop below visits the interior samples one by one and attaches the weight that the pattern $1,4,2,4,\ldots,4,1$ prescribes: an odd index is the midpoint of a parabola and receives weight $4$, an even index is shared between two parabolas and receives weight $2$. The two end points, with weight $1$, are added outside the loop.
def simpson_Mx(N):
# mean-square value computed with Simpson's rule on N points (N odd)
if N % 2 == 0:
raise ValueError("Simpson's Rule requires an odd number of sample points.")
t = np.linspace(ta, tb, N)
y = integrand(t)
h = t[1] - t[0]
total = y[0] + y[-1] # the two end points, weight 1
for i in range(1, N - 1):
if i % 2 == 1:
total = total + 4*y[i] # odd index: parabola midpoint, weight 4
else:
total = total + 2*y[i] # even index: shared end point, weight 2
integral = h/3*total
return integral/T0
The same $N=31$ grid is used, so the two rules consume exactly the same $31$ function values and differ only in the weights they attach to them.
Mx_simp = simpson_Mx(N)
rms_simp = np.sqrt(Mx_simp)
err_simp = abs(Mx_simp - Mx_exact)
print('Simpson mean-square value = %.10f' % Mx_simp)
print('Simpson RMS value = %.10f' % rms_simp)
print('Absolute error = %.6e' % err_simp)
Simpson result
Simpson's Rule also returns the analytical value to round-off,
$$ M_x^{(S)}=9.000000000000, \qquad |M_x^{(S)}-9|\approx7\times10^{-15}. $$
On this grid the two rules are therefore equally accurate, and both are exact for practical purposes. The residual difference between them, at the fifteenth decimal place, reflects only the order in which the floating-point additions happen to be performed.
Both methods approximate the same integral, but they use different approximations of the integrand.
The Trapezoidal Rule uses piecewise linear approximations.
Simpson's Rule uses piecewise quadratic approximations.
The numerical results can now be compared directly with the analytical benchmark.
print('Method Mean-square RMS value Absolute error')
print('-'*72)
print('Analytical %.10f %.10f 0.000000e+00' % (Mx_exact, rms_exact))
print('Trapezoidal %.10f %.10f %.6e' % (Mx_trap, rms_trap, err_trap))
print('Simpson %.10f %.10f %.6e' % (Mx_simp, rms_simp, err_simp))
2.6.1. Summary¶
| Method | Mean-square value | RMS value | Absolute error |
|---|---|---|---|
| Analytical | 9.0000000000 | 3.0000000000 | 0 |
| Trapezoidal | 9.0000000000 | 3.0000000000 | $\approx5.3\times10^{-15}$ |
| Simpson | 9.0000000000 | 3.0000000000 | $\approx7.1\times10^{-15}$ |
On this grid the two rules agree with the analytical result and with each other. Nothing yet distinguishes them, so the more interesting question is what happens when the grid is changed.
Numerical integration is based on a finite set of sample points.
A natural question is:
What happens if we use more sample points?
We repeat the calculation for several odd values of $N$.
Using an odd number of points allows both the Trapezoidal Rule and Simpson's Rule to be evaluated on the same uniform grids.
Ns = [11, 21, 31, 41, 61, 101]
errT = [] # trapezoidal error for each N
errS = [] # Simpson error for each N
for n in Ns:
errT.append(abs(trapz_Mx(n) - Mx_exact))
errS.append(abs(simpson_Mx(n) - Mx_exact))
print(' N subintervals Trapezoidal error Simpson error')
for i in range(len(Ns)):
print('%5d%14d%20.6e%17.6e' % (Ns[i], Ns[i] - 1, errT[i], errS[i]))
The table does not behave as expected. Refining the grid does not reduce the error gradually: the errors are either of order one or at the level of round-off, with nothing in between. Two entries stand out. At $N=11$ both rules are wrong by about $4.1$, which is almost half of the quantity being computed. At $N=21$ the Trapezoidal Rule is exact while Simpson's Rule is wrong by about $1.37$, so on that grid the lower-order method is the accurate one. From $N=31$ onward both are exact.
The error is defined as
$$ E_N=|M_x^{(\mathrm{numerical})}-M_x^{(\mathrm{analytical})}|. $$
Plotting the error on a logarithmic scale makes it easier to observe how rapidly the numerical approximation improves as the grid becomes finer.
Because several of the errors are at the level of round-off, a floor is applied before taking the logarithm; those points mark the limit of double precision, not a measured error.
eps = np.finfo(float).eps # smallest error that floating-point arithmetic allows
plt.figure(figsize = (6, 4))
plt.loglog(Ns, np.maximum(errT, eps), 'o-', label = 'Trapezoidal')
plt.loglog(Ns, np.maximum(errS, eps), 's-', label = 'Simpson')
plt.axhline(eps, color = 'gray', linestyle = '--', label = 'round-off level')
plt.xlabel('Number of sample points N')
plt.ylabel('Absolute error in mean-square value')
plt.title('Numerical integration error')
plt.grid(True, which = 'both')
plt.legend()
plt.show()
Why the error jumps instead of decaying
The explanation lies in the decomposition of Section 2.1. Writing $\omega_0=2\pi/T_0=0.5$ for the fundamental frequency of the interval, the signal is a sum of harmonics $5\omega_0$ and $3\omega_0$, and squaring it produces only the harmonics
$$ 0,\quad 2\omega_0,\quad 6\omega_0,\quad 8\omega_0,\quad 10\omega_0 . $$
The integrand is therefore a trigonometric polynomial containing no frequency above the tenth harmonic, and the average being computed is exactly the constant term of that polynomial.
Now consider what the Trapezoidal Rule does to a single harmonic. Summing $\cos(k\omega_0 t)$ over $n$ equally spaced points covering one full period gives exactly zero, unless $k$ is a multiple of $n$, in which case every sample lands on the same phase of the cosine and the sum returns the amplitude instead of zero. The rule therefore integrates every harmonic exactly, and returns the constant term exactly, provided that none of the harmonics $2,6,8,10$ is a multiple of the number of subintervals $n$. Only $n\in\{2,3,4,5,6,8,10\}$ violates this. With $n=10$, at $N=11$, the tenth harmonic is aliased to a constant and contaminates the answer; every other grid in the table has $n\ge20$ and is exact.
This is a general and useful property: for a smooth periodic integrand sampled over a whole number of periods, the Trapezoidal Rule is not second-order accurate but effectively exact, as soon as the grid resolves the highest frequency present. It is the reason the trapezoidal rule is used for Fourier coefficients and for contour integrals, where more elaborate rules offer nothing.
Simpson's Rule does not inherit the property, and the reason is visible in an algebraic identity. The composite Simpson result on a grid of $n$ subintervals is a combination of two trapezoidal results, one on that grid and one on the coarser grid of $n/2$ subintervals,
$$ S_h=\frac{4T_h-T_{2h}}{3}. $$
The combination is what cancels the leading $h^2$ error term of the trapezoidal rule and raises the order from two to four. But it also drags in the coarse grid: if $T_{2h}$ is contaminated by aliasing, so is $S_h$. At $N=21$ the fine grid has $n=20$ subintervals and is exact, while the coarse grid has $n/2=10$ and aliases the tenth harmonic, so Simpson inherits an error that the Trapezoidal Rule on the same points does not have. The identity is verified numerically below.
print(' N Simpson (4*T_h - T_2h)/3 T_h T_2h')
for n in Ns:
n_sub = n - 1
fine = trapz_Mx(n) # trapezoidal rule with n_sub subintervals
coarse = trapz_Mx(n_sub//2 + 1) # trapezoidal rule with half as many
print('%5d%18.10f%20.10f%16.10f%16.10f' % (n, simpson_Mx(n), (4*fine - coarse)/3, fine, coarse))
The two middle columns agree to every printed digit, and the last two columns show where the $N=21$ failure comes from: the fine grid gives $9$ exactly, the coarse grid gives $4.90$, and their combination $(4\cdot 9-4.90)/3=10.37$ is the value Simpson returns.
The practical lesson is not that Simpson's Rule is unreliable. It is that a higher order of accuracy is a statement about the leading term of an error expansion, and it earns its advantage only when that expansion describes the situation. For a smooth non-periodic integrand, the expansion applies and Simpson is decisively better, as its $O(h^4)$ label promises. For a periodic integrand sampled over whole periods, the expansion is not what governs the error, and the simpler rule is at least as good.
The difference between the two methods can also be understood geometrically.
For the Trapezoidal Rule, the curve is approximated by straight line segments.
For Simpson's Rule, groups of three neighboring points are used to construct a quadratic approximation.
The following figure draws both reconstructions on the same grid, against the exact integrand.
y_dense = integrand(t_dense)
plt.figure(figsize = (6, 6))
# left panel: the trapezoidal rule joins the samples by straight lines
plt.subplot(2, 1, 1)
plt.plot(t_dense, y_dense, 'k', linewidth = 1.5, label = 'exact $x(t)^2$')
plt.plot(t, x2, 'o-', color = 'tab:green', markersize = 4, label = 'piecewise linear')
plt.fill_between(t, x2, alpha = 0.20, color = 'tab:green')
plt.title('Trapezoidal Rule, N = %d' % N)
plt.xlabel('Time t')
plt.ylabel('$x(t)^2$')
plt.grid(True)
plt.legend()
# right panel: Simpson's rule fits one parabola to every pair of subintervals
plt.subplot(2, 1, 2)
plt.plot(t_dense, y_dense, 'k', linewidth = 1.5, label = 'exact $x(t)^2$')
for start in range(0, N - 2, 2):
t_local = t[start:start + 3] # three neighbouring samples
coeff = np.polyfit(t_local, x2[start:start + 3], 2)
t_fine = np.linspace(t_local[0], t_local[-1], 40)
if start == 0:
plt.plot(t_fine, np.polyval(coeff, t_fine), color = 'tab:green', linewidth = 1.5, label = 'piecewise quadratic')
else:
plt.plot(t_fine, np.polyval(coeff, t_fine), color = 'tab:green', linewidth = 1.5)
plt.plot(t, x2, 'ko', markersize = 3)
plt.title("Simpson's Rule, N = %d" % N)
plt.xlabel('Time t')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
The two panels look very different. The piecewise linear reconstruction cuts the top off every peak and misses a visible amount of area, while the parabolas follow the curve closely enough to be hard to distinguish from it. Yet both rules returned the same exact mean square on this grid. The figure is therefore a useful warning: what a quadrature rule reproduces is a number, not a shape. The trapezoidal rule underestimates every peak and overestimates every trough, and over a whole period those opposite errors cancel exactly, which no inspection of the left panel would suggest.
This lab illustrates an important numerical-analysis workflow.
First, an engineering quantity is expressed as an integral:
$$ M_x = \frac{1}{T_0} \int_{-2\pi}^{2\pi}x(t)^2\,dt. $$
Next, the continuous function is sampled on a numerical grid. The Trapezoidal Rule approximates the integral using piecewise linear segments, while Simpson's Rule uses a higher-order polynomial approximation.
The analytical benchmark,
$$ M_x=9, $$
allows the numerical error to be measured directly, and that measurement produced a result worth remembering. On the $N=31$ grid both rules are exact to round-off. On the coarser $N=21$ grid the Trapezoidal Rule is still exact while Simpson's Rule is wrong by more than one unit. The reason is that the integrand is periodic over the integration interval, which makes the Trapezoidal Rule exact as soon as the grid resolves the highest harmonic present, while Simpson's Rule, being a combination of two trapezoidal results on grids of different spacing, is limited by the coarser of the two.
Two conclusions follow, and both generalize beyond this example. A quadrature rule should be chosen against the structure of the integrand, not by its order of accuracy alone. And an analytical benchmark, wherever one is available, is worth the effort of deriving it: without the exact value $M_x=9$, the $N=21$ result would have looked entirely plausible.
Problem: Harmonic forced vibration of a 1-DOF system
Most engineering dynamics reduces to an initial value problem: the governing equation gives the rate of change at the current instant, and the task is to predict the state at a later time.
In this lab we solve the forced vibration of a single-degree-of-freedom oscillator with the simplest possible time integrator, the explicit Euler method, and examine how it behaves. Because the system is linear, its exact solution can be written down, so every numerical result can be measured against the truth rather than against another approximation.
The lab is organized around one number, the step size $\Delta t$, and the two questions that can be asked about it. Is there a value above which the method fails altogether? There is, and Section 3.7 derives it: $\Delta t=0.1$ s for this oscillator, below which the computed motion converges to the true one and above which it grows without bound while the true motion decays. And once the step is below that limit, how does the computed motion approach the true one as the step is refined? Section 3.4 answers that from the formula itself.
Consider a mass $m$ supported by a spring of stiffness $k$ and a viscous damper of coefficient $c$, driven by a harmonic force:
$$ m\ddot{x}+c\dot{x}+kx=F_0\sin(\omega t), \qquad x(0)=x_0,\quad \dot{x}(0)=v_0. $$
From a second-order equation to state-space form
Numerical integrators are written for first-order systems, so the equation of motion must be rewritten before it can be marched. Introducing the state vector
$$ \mathbf{y}=\begin{bmatrix}x\\v\end{bmatrix}, \qquad v=\dot{x}, $$
the second-order equation becomes two coupled first-order equations,
$$ \dot{\mathbf{y}} = \begin{bmatrix} v\\ \dfrac{F_0\sin\omega t-cv-kx}{m} \end{bmatrix}. $$
The first row is the definition of velocity and the second is Newton's second law solved for acceleration. The conversion is general: an $n$-th order equation becomes a system of $n$ first-order equations, which is why solver libraries provide first-order integrators only.
The physical parameters are chosen so that the response contains both a decaying transient and a steady oscillation.
| Quantity | Symbol | Value |
|---|---|---|
| Mass | $m$ | $1.0$ kg |
| Damping coefficient | $c$ | $0.4$ N·s/m |
| Stiffness | $k$ | $4.0$ N/m |
| Force amplitude | $F_0$ | $1.0$ N |
| Forcing frequency | $\omega$ | $1.2$ rad/s |
| Initial displacement | $x_0$ | $0.5$ m |
| Initial velocity | $v_0$ | $0$ m/s |
| Final time | $T$ | $20$ s |
These values give a natural frequency $\omega_n=\sqrt{k/m}=2$ rad/s and a damping ratio $\zeta=c/(2\sqrt{km})=0.1$, so the system is lightly damped. The forcing frequency $\omega=1.2$ rad/s lies below resonance. The response therefore superposes a transient at the damped natural frequency $\omega_d=\omega_n\sqrt{1-\zeta^2}\approx1.99$ rad/s on a steady oscillation at $1.2$ rad/s, and the integrator must capture both.
The symbols are those of the equation of motion. The derived quantities $\omega_n$, $\zeta$, $\omega_d$ are used by the exact solution in Section 3.5 and by the stability analysis in Section 3.7.
# --- Oscillator ---
m = 1.0 # mass [kg]
c = 0.4 # damping [N s/m]
k = 4.0 # stiffness [N/m]
F0 = 1.0 # force amplitude [N]
w = 1.2 # forcing frequency [rad/s]
x0 = 0.5 # initial displacement [m]
v0 = 0.0 # initial velocity [m/s]
T = 20.0 # final time [s]
wn = np.sqrt(k/m) # natural frequency
zeta = c/(2*np.sqrt(k*m)) # damping ratio
wd = wn*np.sqrt(1 - zeta**2) # damped frequency
print('wn = %.3f rad/s, zeta = %.3f, wd = %.3f rad/s' % (wn, zeta, wd))
The integrator evaluates one thing at every step: the acceleration, which Newton's second law gives from the current time, displacement and velocity,
$$ a(t,x,v)=\frac{F_0\sin\omega t-cv-kx}{m}. $$
Together with the definition $\dot x=v$ this is the right-hand side of the state-space system. Writing it once, separately from the integrator, is what makes the marching loop independent of the problem: replacing this one function turns the same loop into a solver for a different oscillator.
def acceleration(t, x, v):
# Newton's second law solved for the acceleration
return (F0*np.sin(w*t) - c*v - k*x)/m
Expand the solution about the current time and keep two terms,
$$ \mathbf{y}(t+\Delta t) = \mathbf{y}(t)+\Delta t\,\dot{\mathbf{y}}(t) +\frac{\Delta t^2}{2}\ddot{\mathbf{y}}(t)+\cdots, $$
then replace $\dot{\mathbf{y}}$ by the right-hand side of the differential equation and discard everything beyond the linear term:
$$ \mathbf{y}_{n+1}=\mathbf{y}_n+\Delta t\,\mathbf{f}(t_n,\mathbf{y}_n). $$
Geometrically the method evaluates the slope at the current point and follows the tangent line for one step. Where the true solution curves away from that line the step lands slightly off, and the next step starts from the wrong place.
The first discarded term carries $\Delta t^2$, so a single step is wrong by $O(\Delta t^2)$. What matters in practice is the error at the end of the simulation, and these are not the same. Reaching a fixed final time $T$ requires $T/\Delta t$ steps, so halving the step halves each individual mistake but doubles how many are made:
$$ \frac{T}{\Delta t}\times O(\Delta t^2)=O(\Delta t). $$
One order is lost in the accumulation, which is why the Euler method is called first-order accurate. The two runs of Section 3.6 are both consistent with it.
One function marches the whole simulation: it builds the time grid, stores the displacement and the velocity in two arrays with one entry per time level, sets the initial state, and applies the two update formulas in a loop, first for $x$ and then for $v$. Both arrays are returned, although the rest of the lab uses only the displacement.
def euler(dt):
t = np.arange(0, T + dt/2, dt)
nt = len(t)
x = np.zeros(nt) # x[n] = displacement at t[n]
v = np.zeros(nt) # v[n] = velocity at t[n]
x[0] = x0
v[0] = v0
for n in range(nt - 1):
a = acceleration(t[n], x[n], v[n]) # slope at the current point ...
x[n+1] = x[n] + dt*v[n] # ... followed for one step
v[n+1] = v[n] + dt*a
return t, x, v
Measuring an error requires something to measure against. For this linear system the exact solution is available, which makes the comparison unambiguous: every discrepancy below is the fault of the integrator and nothing else.
The response is the sum of a particular solution at the forcing frequency and a homogeneous solution that decays,
$$ x(t) = \underbrace{e^{-\zeta\omega_n t}\left(A\cos\omega_d t+B\sin\omega_d t\right)}_{\text{transient}} + \underbrace{X\sin(\omega t-\varphi)}_{\text{steady state}}, $$
with the steady-state amplitude and phase
$$ X=\frac{F_0}{\sqrt{(k-m\omega^2)^2+(c\omega)^2}}, \qquad \varphi=\arctan\frac{c\omega}{k-m\omega^2}, $$
and the constants $A$ and $B$ fixed by the initial conditions:
$$ A=x_0+X\sin\varphi, \qquad B=\frac{v_0-X\omega\cos\varphi+\zeta\omega_n A}{\omega_d}. $$
The transient decays by the factor $e^{-\zeta\omega_n t}$, so by $t=20$ s it should have fallen to $e^{-4}\approx1.8\%$ of its initial amplitude. That number is the benchmark for Section 3.6: a computed transient that has not decayed by then is wrong.
X = F0/np.sqrt((k - m*w**2)**2 + (c*w)**2)
phi = np.arctan2(c*w, k - m*w**2)
A = x0 + X*np.sin(phi)
B = (v0 - X*w*np.cos(phi) + zeta*wn*A)/wd
def exact(t):
return np.exp(-zeta*wn*t)*(A*np.cos(wd*t) + B*np.sin(wd*t)) + X*np.sin(w*t - phi)
print('X = %.4f m, phi = %.4f rad, A = %.4f, B = %.4f' % (X, phi, A, B))
With the method and the reference in place, the first question is simply whether the computed motion is right. Two step sizes are tried, one small and one on the other side of the limit that Section 3.7 will derive: $\Delta t=0.01$ s and $\Delta t=0.12$ s.
td = np.linspace(0, T, 1000)
t1, x1, v1 = euler(0.01)
t2, x2, v2 = euler(0.12)
plt.figure(figsize = (6, 4))
plt.plot(td, exact(td), 'k', linewidth = 5, alpha = 0.3, label = 'exact')
plt.plot(t1, x1, 'b--', label = 'Euler, dt = 0.01')
plt.plot(t2, x2, 'r--', label = 'Euler, dt = 0.12')
plt.xlim([0, 20])
plt.xlabel('t [s]')
plt.ylabel('x [m]')
plt.legend()
plt.show()
print('exact : max |x| = %.3f m' % abs(exact(td)).max())
print('dt = 0.01 : max |x| = %.3f m, RMS error = %.3e m' % (abs(x1).max(), np.sqrt(np.mean((x1 - exact(t1))**2))))
print('dt = 0.12 : max |x| = %.3f m, RMS error = %.3e m' % (abs(x2).max(), np.sqrt(np.mean((x2 - exact(t2))**2))))
What the figure shows
At $\Delta t=0.01$ s the Euler curve follows the exact solution closely enough that the two are hard to tell apart, parting from it only slightly at the largest peaks. The root-mean-square error is $1.1\times10^{-2}$ m, about two percent of the peak amplitude. The method works, and by the first-order argument of Section 3.4 the remaining error would halve again with each halving of the step.
At $\Delta t=0.12$ s the same method applied to the same equation does something the true motion never does: it grows. The exact transient decays by the factor $e^{-4}\approx0.018$ over the run, so that by $t=20$ s only the steady forced oscillation of amplitude $0.38$ m remains. The computed oscillation instead swells on every cycle and reaches $1.5$ m at the end of the run, four times the true value and still growing. Nothing about the forcing, the initial condition or the physical parameters has changed between the two curves. Only the step size has.
The mechanism is specific. One Euler step multiplies the free-vibration part of the state by the factor $\lvert 1+\lambda\Delta t\rvert$, where $\lambda$ is an eigenvalue of the system; the exact solution multiplies it by $e^{-\zeta\omega_n\Delta t}$ over the same interval. At $\Delta t=0.12$ s the Euler factor is $1.0048$, larger than one, so each step adds a little energy that the damper should have removed, and over the $167$ steps of the run $1.0048^{167}\approx2.2$. The computed transient is amplified where the true one decays, and that is what the figure shows. Following the tangent line always carries the state slightly outside the true trajectory; whether the error shrinks or compounds depends only on whether the factor is below or above one, and that is the question the next section answers.
Set the forcing aside for a moment. Errors do not obey the forced equation; they obey the homogeneous one,
$$ \dot{\mathbf{y}}=A\mathbf{y}, \qquad A=\begin{bmatrix}0&1\\-k/m&-c/m\end{bmatrix} =\begin{bmatrix}0&1\\-4&-0.4\end{bmatrix}. $$
Because the eigenvectors of $A$ decouple the system into independent scalar equations, the scalar analysis applies to each eigenvalue $\lambda$ separately. One explicit Euler step multiplies each component by the amplification factor $1+\lambda\Delta t$, so the method is stable only if
$$ \left\lvert 1+\lambda\Delta t\right\rvert<1 \qquad \text{for every eigenvalue } \lambda . $$
Give the product a name, $z=\lambda\Delta t$. It is a complex number, one per eigenvalue, and the condition reads $\lvert 1+z\rvert<1$: the distance from $z$ to the point $-1$ must be less than one. In the complex plane this says that $z$ must lie inside the unit disk centred at $-1$, which is called the stability region of the method. A step size does not move the eigenvalues; it scales $z$ along the ray from the origin through $\lambda$, so a small $\Delta t$ keeps $z$ near the origin, which is on the boundary of the disk, and a larger $\Delta t$ pushes it outward along that ray until it leaves the disk. Writing $\lambda=-\zeta\omega_n\pm i\omega_d$ and expanding the inequality gives
$$ \Delta t<\frac{2\zeta\omega_n}{\lvert\lambda\rvert^2}=\frac{2\times0.2}{4}=0.1\ \text{s}. $$
The two runs of Section 3.6 sit on either side of this number, $0.01$ s well inside and $0.12$ s just outside, which is exactly why one converged and the other grew. The cells below compute $z$ for both and draw them on the disk.
Two features of this condition deserve attention. It involves only the properties of the system, not the accuracy the user wants. And the forcing amplitude $F_0$ does not appear in it at all, because stability concerns the growth of errors rather than the size of the response.
A_matrix = np.array([[0, 1], [-k/m, -c/m]]) # the system written as y' = A y
lam = np.linalg.eigvals(A_matrix)
print('eigenvalues lam =', np.round(lam, 3))
print('stability limit dt < %.3f s' % (-2*lam[0].real/abs(lam[0])**2))
print()
for dt in [0.01, 0.12]:
z = lam[0]*dt # the upper eigenvalue; the lower one is its mirror image
print('dt = %.2f : z = lam dt = %.3f %+.3fi |1 + z| = %.4f' % (dt, z.real, z.imag, abs(1 + z)))
# the stability region: the unit disk centred at -1, |1 + z| < 1
th = np.linspace(0, 2*np.pi, 200)
z1 = lam*0.01
z2 = lam*0.12
plt.figure(figsize = (6, 8))
# left: the whole disk
plt.subplot(2, 1, 1)
plt.plot(-1 + np.cos(th), np.sin(th), 'k', label = '|1 + z| = 1')
plt.plot(-1, 0, 'k+', markersize = 12)
plt.plot(z1.real, z1.imag, 'bo', label = 'z = lam dt, dt = 0.01')
plt.plot(z2.real, z2.imag, 'ro', label = 'z = lam dt, dt = 0.12')
plt.plot([-1, z2[0].real], [0, z2[0].imag], 'r:') # distance |1 + z| from the centre
plt.text(-1, -0.12, '-1', ha = 'center')
plt.axis('equal')
plt.xlabel('Re z')
plt.ylabel('Im z')
plt.title('Stability region of explicit Euler')
plt.legend(loc = 'lower left')
# right: zoom on the origin, where both z sit
plt.subplot(2, 1, 2)
plt.plot(-1 + np.cos(th), np.sin(th), 'k', label = '|1 + z| = 1')
plt.plot(0, 0, 'k+', markersize = 12)
plt.plot(z1.real, z1.imag, 'bo')
plt.plot(z2.real, z2.imag, 'ro')
plt.text(z1[0].real + 0.01, z1[0].imag, 'z = %.3f + %.3fi' % (z1[0].real, z1[0].imag), color = 'b')
plt.text(z2[0].real + 0.01, z2[0].imag, 'z = %.3f + %.3fi' % (z2[0].real, z2[0].imag), color = 'r')
plt.xlim(-0.12, 0.12)
plt.ylim(-0.3, 0.3)
plt.xlabel('Re z')
plt.ylabel('Im z')
plt.title('Near the origin: inside (blue), outside (red)')
plt.legend(loc = 'lower left')
plt.tight_layout()
plt.show()
Stability result
The printed factors settle the question: $\lvert 1+z\rvert=0.9982$ at $\Delta t=0.01$ s and $1.0048$ at $\Delta t=0.12$ s, one below and one above one. The figure says the same thing geometrically. The left panel shows the whole stability region, the unit disk centred at $-1$, and the dotted line is the distance $\lvert 1+z\rvert$ from that centre to the red point; it is just longer than the radius. Both pairs of points sit close to the origin, which lies on the boundary of the disk, so the right panel zooms in. There the blue pair, $z=-0.002\pm0.020\,i$, is just inside the arc, and the red pair, $z=-0.024\pm0.239\,i$, is just outside it. The two points lie on the same ray from the origin, because $z=\lambda\Delta t$ and only $\Delta t$ changed; the larger step has simply carried $z$ across the boundary.
Notice how close to the boundary even the stable case is. The eigenvalues are almost purely imaginary, so the ray from the origin runs almost tangent to the circle, $\lvert 1+z\rvert$ is only slightly below one for any permitted step, and the numerical damping per step is very weak. An explicit Euler step applied to a lightly damped oscillator is barely contractive inside the limit, and above $\Delta t=0.1$ s it is not contractive at all. The limit is a threshold, not a guideline: at $0.12$ s the growth factor per step is only $1.0048$, less than half a percent, and it was still enough to multiply the transient by $2.2$ over a run in which it should have fallen to $0.018$. A larger step would only make the same failure arrive sooner.
This lab used one method in order to separate three things that are easy to confuse.
Stability is whether errors decay or grow as the marching proceeds. It is governed by the eigenvalues of the system and the step size, through the condition $\lvert 1+\lambda\Delta t\rvert<1$, and it does not depend on the accuracy the user wants or on the forcing. For this oscillator the condition is $\Delta t<0.1$ s. Below it the computed motion converges to the true one; above it the computed motion grows while the true one decays, and no amount of further marching improves anything.
Accuracy is how closely a convergent numerical solution follows the true one. It is governed by the order of the method and improves predictably as the step is refined, in proportion to $\Delta t$ for the explicit Euler method, as the Taylor argument of Section 3.4 showed.
Cost is the number of right-hand-side evaluations, one per step here, and it is the price paid for the accuracy obtained.
For this lightly damped oscillator the three came together in one conclusion. The eigenvalues lie close to the boundary of the Euler stability region, so the method is simultaneously barely stable and slowly convergent: at $\Delta t=0.12$ s the solution diverges, and at $\Delta t=0.01$ s it reproduces the motion to within two percent, so that two thousand steps buy barely two correct digits. Both limitations trace back to the same source, a single slope evaluated at the beginning of the interval, and both are what the higher-order and implicit methods of a fuller treatment are built to address.
Problem: Transverse vibration of a stretched string
The wave equation is the second of the two classical partial differential equations of engineering, and it behaves in almost the opposite way to the first. Diffusion smooths a profile out and forgets its initial shape; a wave carries the shape along at a finite speed and, in the absence of damping, never loses it.
In this lab we solve the vibration of a string fixed at both ends with the standard explicit scheme, central differences in both time and space, and examine the one condition that governs whether the computation succeeds at all.
The lab is built around a result that has no counterpart in the heat equation:
- the scheme is stable only if the Courant number $C=c\,\Delta t/\Delta x$ satisfies $C\le1$,
- the accuracy improves as $C$ approaches that limit,
- at exactly $C=1$ the scheme reproduces the exact solution to round-off.
Taking a smaller time step therefore makes this computation slower and less accurate, which is worth understanding before it is met by accident.
A string of length $L$ under tension, fixed at both ends, has a transverse displacement $u(x,t)$ governed by
$$ \frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}, \qquad 0<x<L,\quad t>0, $$
where $c=\sqrt{T_s/\rho}$ is the wave speed, $T_s$ the tension and $\rho$ the mass per unit length. The ends are held fixed,
$$ u(0,t)=u(L,t)=0, $$
and because the equation is second order in time, two initial conditions are required:
$$ u(x,0)=f(x), \qquad \frac{\partial u}{\partial t}(x,0)=g(x). $$
The string is released from rest in this lab, so $g(x)=0$.
What makes this different from diffusion. The heat equation is first order in time and its solutions decay: every Fourier mode is multiplied by $e^{-\alpha k^2 t}$, so high wavenumbers vanish quickly and the solution becomes smoother as time passes. The wave equation is second order in time and its solutions oscillate: every mode is multiplied by $\cos(kct)$ or $\sin(kct)$, so nothing decays, nothing smooths, and an error introduced at any point is carried forward forever. A numerical scheme for the wave equation therefore has no diffusion to hide its mistakes, which is why the stability condition below is so sharp.
The numbers are chosen so that the exact solution is a single standing wave and the period is convenient.
| Quantity | Symbol | Value |
|---|---|---|
| Wave speed | $c$ | $1.0$ m/s |
| String length | $L$ | $1.0$ m |
| Final time | $T$ | $5.0$ s |
| Grid points | $n_x$ | $41$ |
| Initial shape | $f(x)$ | $\sin(\pi x)$ |
| Initial velocity | $g(x)$ | $0$ |
With $c=L=1$ the fundamental period is $2L/c=2$ s, so the final time covers two and a half periods. The grid of $41$ points gives $\Delta x=0.025$ m; the time step is not chosen directly but through the Courant number, $\Delta t=C\,\Delta x/c$.
The time step is absent from the list below. It is set later by the Courant number, which is the quantity this lab is about.
c = 1.0 # wave speed [m/s]
L = 1.0 # string length [m]
T = 5.0 # final time [s]
nx = 41 # grid points
print('dx = %.4f m, fundamental period = %.1f s' % (L/(nx - 1), 2*L/c))
Place a grid $x_i=i\,\Delta x$ and time levels $t_n=n\,\Delta t$, and write $u_i^n\approx u(x_i,t_n)$. Both derivatives in the equation are second derivatives, so both are replaced by the same three-point central formula:
$$ \frac{\partial^2 u}{\partial t^2}\bigg|_i^n \approx \frac{u_i^{n+1}-2u_i^{n}+u_i^{n-1}}{\Delta t^2}, \qquad \frac{\partial^2 u}{\partial x^2}\bigg|_i^n \approx \frac{u_{i+1}^{n}-2u_i^{n}+u_{i-1}^{n}}{\Delta x^2}. $$
Substituting both into the wave equation and solving for the only unknown, the value at the new time level,
$$ u_i^{n+1} = 2u_i^{n}-u_i^{n-1} +C^2\left(u_{i+1}^{n}-2u_i^{n}+u_{i-1}^{n}\right), \qquad C=\frac{c\,\Delta t}{\Delta x}. $$
The dimensionless group $C$ is the Courant number. It is the ratio of the distance a physical wave travels in one time step, $c\,\Delta t$, to the width of one grid cell, $\Delta x$, and it is the only parameter the scheme depends on. Everything in this lab is a statement about $C$.
Two structural remarks are worth making before the code.
The scheme is explicit: the new value is computed directly from known values, with no system to solve, exactly like the FTCS scheme for the heat equation. But it uses three time levels rather than two, because the time derivative is second order. That means the marching loop must carry two arrays of history, and it means the first step cannot be taken by the formula above.
The stencil reaches only one cell to each side. In one time step, information in the numerical solution can therefore travel at most one grid spacing, whatever the physics demands. That observation is all that is needed for the stability condition of Section 4.8.
At $n=0$ the update formula requires $u_i^{-1}$, a value one step before the initial condition, which does not exist. The second initial condition supplies it. Approximating $u_t(x,0)=g(x)$ with a central difference in time,
$$ \frac{u_i^{1}-u_i^{-1}}{2\Delta t}=g_i \qquad\Longrightarrow\qquad u_i^{-1}=u_i^{1}-2\Delta t\,g_i , $$
and substituting this into the update formula written at $n=0$ eliminates the fictitious level and gives a special starting formula,
$$ u_i^{1} = u_i^{0}+\Delta t\,g_i +\frac{C^2}{2}\left(u_{i+1}^{0}-2u_i^{0}+u_{i-1}^{0}\right). $$
With the string released from rest, $g=0$, the middle term drops out.
Remark. It is tempting to start instead with $u_i^1=u_i^0$, on the grounds that the string has zero initial velocity and therefore does not move during the first instant. That is a first-order approximation, and using it degrades the entire computation to first order no matter how accurate the scheme is afterwards: refining the mesh by a factor of two then reduces the error by two rather than by four. A single sloppy step at the start is enough to destroy the order of accuracy of every step that follows.
The whole method is one function. It builds the grid, stores the entire history in one array U with one row per time level and one column per grid point, so that U[n, i] is $u_i^n$ exactly as in the formulas, takes the special first step, then marches the three-level formula with one loop over time levels and one over interior points. The two end columns are never written after the initial condition, which is how the fixed boundaries are enforced. The full history is returned so that the sections below can plot any instant.
def wave(C, nx):
x = np.linspace(0, L, nx)
dx = x[1] - x[0]
dt = C*dx/c
nt = int(round(T/dt)) # number of time steps
t = np.arange(nt + 1)*dt
U = np.zeros((nt + 1, nx)) # U[n, i] = u at x_i and t_n
U[0, :] = np.sin(np.pi*x) # initial shape, released from rest
# special first step (g = 0)
for i in range(1, nx - 1):
U[1, i] = U[0, i] + 0.5*C**2*(U[0, i+1] - 2*U[0, i] + U[0, i-1])
# three-level scheme
for n in range(1, nt):
for i in range(1, nx - 1):
U[n+1, i] = 2*U[n, i] - U[n-1, i] + C**2*(U[n, i+1] - 2*U[n, i] + U[n, i-1])
return x, t, U
The initial shape $f(x)=\sin(\pi x)$ is the fundamental mode of the fixed string, and a single mode released from rest simply oscillates in place:
$$ u(x,t)=\sin(\pi x)\cos(\pi c\,t). $$
Substituting confirms it: differentiating twice in time gives $-(\pi c)^2 u$ and twice in space gives $-\pi^2 u$, so $u_{tt}=c^2u_{xx}$ holds exactly. The boundary conditions hold because $\sin(0)=\sin(\pi)=0$, and the initial conditions because $\cos(0)=1$ and the time derivative carries a $\sin(\pi c t)$ that vanishes at $t=0$.
This is a standing wave: the shape never changes, only its amplitude, which swings between $+1$ and $-1$ with period $2L/c=2$ s. Having it in closed form means every discrepancy below is the fault of the scheme and nothing else.
def exact(x, t):
return np.sin(np.pi*x)*np.cos(np.pi*c*t)
Run the scheme at $C=0.5$, safely below the limit derived in the next section, and compare it with the exact standing wave.
x, t, U = wave(0.5, nx)
dt = t[1] - t[0]
plt.figure(figsize = (6, 4))
for ti in [0.0, 0.3, 0.6, 0.9]:
n = int(round(ti/dt)) # time level closest to t = ti
if ti == 0.0:
plt.plot(x, exact(x, t[n]), 'k', label = 'exact', lw = 3, alpha = 0.3)
plt.plot(x, U[n, :], 'r--', label = 'CTCS')
else:
plt.plot(x, exact(x, t[n]), 'k', lw = 3, alpha = 0.3)
plt.plot(x, U[n, :], 'r--')
plt.xlabel('x [m]')
plt.ylabel('u [m]')
plt.legend()
plt.show()
print('C = 0.5 : dt = %.4f s, %d steps' % (dt, len(t) - 1))
The two sets of curves lie on top of one another: the scheme reproduces the standing wave, its shape, its amplitude and its rhythm. At this resolution the error is not visible, so it has to be measured.
The curves above are four slices of one surface. Because the whole history is stored as U[n, i], the solution can also be drawn over the $(x,t)$ plane as a surface, here for the first period, $0\le t\le 2$ s, with the same four instants marked on it in red. The surface makes the standing-wave character visible at a glance: every point of the string moves up and down in step, the nodes at the two ends never move, and the amplitude swings between $+1$ and $-1$ with the period of $2$ s.
x, t, U = wave(0.5, nx)
dt = t[1] - t[0]
n_period = int(round(2*L/c/dt)) # time levels in one period (2 s)
X, Tg = np.meshgrid(x, t[:n_period + 1]) # X[n, i] = x_i, Tg[n, i] = t_n, U[n, i] = u(x_i, t_n)
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Tg, U[:n_period + 1, :], cmap = 'viridis', rstride = 4, cstride = 1, linewidth = 0, alpha = 0.5)
for ti in [0.0, 0.3, 0.6, 0.9]:
n = int(round(ti/dt)) # time level closest to t = ti
ax.plot(x, ti*np.ones(nx), U[n, :], linewidth = 2, color = 'r')
ax.set_xlabel('x [m]', fontsize = 12)
ax.set_ylabel('t [s]', fontsize = 12)
ax.set_zlabel('u [m]', fontsize = 12)
ax.view_init(elev = 30, azim = -45)
plt.show()
plt.figure(figsize = (6, 4))
for C in [0.5, 1.0]:
x, t, U = wave(C, nx)
err = np.zeros(len(t))
for n in range(len(t)):
err[n] = np.abs(U[n, :] - exact(x, t[n])).max()
plt.semilogy(t, np.maximum(err, 1e-17), label = 'C = %.1f' % C)
plt.xlabel('t [s]')
plt.ylabel('max |error| [m]')
plt.legend()
plt.show()
What the figure shows
At $C=0.5$ the error grows steadily with time, reaching about $3\times10^{-3}$ after two and a half periods. The growth is linear rather than explosive, which is the signature of a phase error: the numerical wave oscillates at a very slightly wrong frequency, so it drifts out of step with the exact one a little more on every cycle. Section 4.9 computes that frequency.
At $C=1$ the error sits at $10^{-16}$ for the whole simulation. That is round-off, and it does not grow. The scheme has not approximated the solution at all; it has reproduced it. Section 4.9 explains why.
Substitute a single Fourier mode $u_i^n=A^n e^{\,\mathrm{i}kx_i}$ into the update formula, where $A$ is the factor by which the mode is multiplied each step. Using $e^{\,\mathrm{i}k\Delta x}-2+e^{-\mathrm{i}k\Delta x}=-4\sin^2(k\Delta x/2)$,
$$ A-2+A^{-1}=-4C^2\sin^2\!\left(\frac{k\Delta x}{2}\right), $$
which rearranges into a quadratic for the amplification factor,
$$ A^2-2\left(1-2s^2\right)A+1=0, \qquad s=C\sin\!\left(\frac{k\Delta x}{2}\right). $$
The product of the two roots is the constant term, $1$. That single fact decides everything. If the roots are complex they must then be a conjugate pair of modulus exactly one, so the mode neither grows nor decays, which is precisely what a wave should do. If the roots are real their product is still one, so one of them has modulus greater than one and that mode grows without bound.
The roots are complex when $\left|1-2s^2\right|<1$, that is when $s<1$. At $s=1$ exactly they merge into the double root $A=-1$, still of modulus one, and beyond it they separate along the real axis. Since $\sin(k\Delta x/2)$ reaches $1$ for the shortest wave the grid can carry, the condition must hold for $s=C$, and the requirement on the scheme is
$$ C=\frac{c\,\Delta t}{\Delta x}\le 1 $$
This is the CFL condition, after Courant, Friedrichs and Lewy.
The same condition, physically. The stencil of the scheme reaches one cell to each side per step, so the numerical solution at a point can only be influenced by data within a distance $\Delta x$ per step. The true solution at that point is influenced by data within $c\,\Delta t$, since that is how far the wave travels. If $c\,\Delta t>\Delta x$, the exact answer depends on information the scheme never looks at, and no amount of arithmetic can recover it. The numerical domain of dependence must contain the physical one. This argument is older than the Fourier analysis and applies to schemes for which the analysis is intractable.
Notice what the condition is not. For the explicit heat equation the limit is $\alpha\Delta t/\Delta x^2\le1/2$, so halving the mesh forces a fourfold reduction in the time step. Here the restriction is linear in $\Delta x$: halving the mesh halves the permitted time step, no worse. Wave problems are much cheaper to march explicitly than diffusion problems for this reason.
The factor is evaluated for the shortest wave the grid can carry, the one that alternates sign from point to point, for which $\sin(k\Delta x/2)=1$. If any wave is amplified it is that one.
for C in [0.5, 1.0, 1.02]:
s = C # the worst case, sin(k dx / 2) = 1
p = 2*(1 - 2*s**2) # A^2 - p A + 1 = 0
if p**2 < 4:
print('C = %4.2f : complex roots, |A| = 1.0000 -> neutral' % C)
elif p**2 == 4:
print('C = %4.2f : double root, |A| = 1.0000 -> marginal' % C)
else:
amp = abs((p - np.sqrt(p**2 - 4))/2)
print('C = %4.2f : real roots, |A| = %.4f -> growth' % (C, amp))
The running maximum is followed rather than the instantaneous one, because the true amplitude passes through zero twice per period, when the string is momentarily flat, and those dips would obscure the comparison.
plt.figure(figsize = (6, 4))
for C in [0.5, 1.0, 1.02]:
x, t, U = wave(C, nx)
peak = np.zeros(len(t)) # largest |u| seen up to each time level
peak[0] = np.abs(U[0, :]).max()
for n in range(1, len(t)):
peak[n] = max(peak[n-1], np.abs(U[n, :]).max())
plt.semilogy(t, peak, label = 'C = %.2f' % C)
plt.xlabel('t [s]')
plt.ylabel('largest |u| so far [m]')
plt.legend()
plt.show()
x, t, U = wave(1.02, nx)
print('C = 1.02 : max |u| = %.2e m at t = %.2f s' % (np.abs(U[-1, :]).max(), t[-1]))
Stability result
At $C=0.5$ and $C=1$ the two curves lie on top of one another at exactly $1$ for the whole simulation, as they should: the string is undamped and the true solution neither grows nor decays. At $C=1.02$, only two percent above the limit, the displacement reaches $10^{16}$ m.
The detail worth noticing is where that came from. The initial condition is a single smooth mode with no short-wavelength content whatsoever. The unstable mode, the sawtooth that alternates sign from one grid point to the next, is present only at the level of round-off, around $10^{-16}$. But it is multiplied by $1.49$ on every step, and the simulation takes about $200$ steps: $1.49^{200}$ is roughly $10^{34}$. An error the size of the last bit of a floating-point number is enough, and this is why a stability violation cannot be tolerated for even a short run. The curve on the plot shows exactly that history, flat while the instability is still below the physical solution, then a straight line on the logarithmic axis once it emerges.
The same substitution that produced the stability condition also gives the frequency at which the numerical solution oscillates. Writing the complex roots as $A=e^{\pm\mathrm{i}\omega\Delta t}$ and matching real parts, $\cos(\omega\Delta t)=1-2s^2$, which simplifies to
$$ \sin\!\left(\frac{\omega\,\Delta t}{2}\right)=C\sin\!\left(\frac{k\Delta x}{2}\right) \qquad\Longrightarrow\qquad \omega_{\text{num}}=\frac{2}{\Delta t}\arcsin\!\left[C\sin\!\left(\frac{k\Delta x}{2}\right)\right]. $$
The exact wave equation has $\omega=ck$. The numerical scheme does not, and the mismatch is called numerical dispersion: different wavenumbers travel at slightly different speeds on the grid, whereas in the true equation they all travel at $c$. For a single mode, as here, the only visible consequence is the slow phase drift seen in Section 4.7.
Now set $C=1$ in that formula. Then $\Delta x=c\,\Delta t$, and
$$ \omega_{\text{num}}=\frac{2}{\Delta t}\arcsin\!\left[\sin\!\left(\frac{k\Delta x}{2}\right)\right] =\frac{k\Delta x}{\Delta t}=kc , $$
which is exact, for every wavenumber at once. The reason is visible in the update formula itself. Putting $C=1$ makes the $u_i^n$ terms cancel completely:
$$ u_i^{n+1}=u_{i+1}^{n}+u_{i-1}^{n}-u_i^{n-1}. $$
This is d'Alembert's solution written on the grid. At $C=1$ the grid diagonals coincide with the characteristics $x\pm ct=\text{constant}$, so the scheme does not approximate the propagation at all: it simply copies values along the characteristics, which is what the exact solution does. The consequence is the one announced at the start of the lab, and it inverts the usual intuition.
Cs = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1.0]
err = []
for C in Cs:
x, t, U = wave(C, nx)
worst = 0.0
for n in range(len(t)):
worst = max(worst, np.abs(U[n, :] - exact(x, t[n])).max())
err.append(worst)
plt.figure(figsize = (6, 4))
plt.semilogy(Cs, np.maximum(err, 1e-17), 'ro-')
plt.xlabel('Courant number C')
plt.ylabel('max |error| [m]')
plt.show()
for i in range(len(Cs)):
steps = round(T/(Cs[i]*L/(nx - 1)/c))
print('C = %4.2f steps = %5d max error = %.3e' % (Cs[i], steps, err[i]))
Reading the table
The error falls monotonically as $C$ rises, and at $C=1$ it drops by fourteen orders of magnitude to round-off. A smaller time step costs more work and delivers a worse answer: at $C=0.2$ the computation takes five times as many steps as at $C=1$ and is wrong by $3.5\times10^{-3}$ instead of being right.
This is the opposite of the behaviour of the explicit heat equation, where the stability limit $r\le1/2$ is an obstacle to be kept away from and refining the time step always helps. For the wave equation the stability limit is a target. In practice one runs at $C$ just below $1$, both for accuracy and because it minimizes the number of steps, and a safety margin of a few percent is kept only because $c$ may vary over the domain in a real problem.
The one thing that does behave conventionally is mesh refinement at a fixed Courant number, checked next. Because $\Delta t$ and $\Delta x$ are then refined together, both central differences improve at the same rate and the scheme shows its second-order accuracy.
prev = 0.0
for n in [21, 41, 81, 161]:
x, t, U = wave(0.5, n)
worst = 0.0
for j in range(len(t)):
worst = max(worst, np.abs(U[j, :] - exact(x, t[j])).max())
if prev == 0.0:
print('nx = %4d dx = %.5f error = %.3e' % (n, x[1] - x[0], worst))
else:
print('nx = %4d dx = %.5f error = %.3e ratio = %.2f' % (n, x[1] - x[0], worst, prev/worst)) # halving dx should quarter the error
prev = worst
The ratio is $4$ at every refinement, confirming that the scheme is second-order accurate in $\Delta x$, and therefore in $\Delta t$ as well since the two are tied together by $C$. Halving the mesh costs four times the work, two for the extra grid points and two for the extra time steps, and buys a factor of four in accuracy.
Three properties of this scheme are worth carrying away, and only the first is shared with the diffusion problems met earlier.
It is second-order accurate in both variables. Both derivatives were replaced by the same three-point central formula, each with an $O(h^2)$ truncation error, and the mesh refinement study confirms the rate. Getting that rate also required the first time step to be taken with the special second-order formula of Section 4.4; the obvious shortcut would have silently cost an order.
It is conditionally stable, with a linear condition. The CFL condition $c\,\Delta t\le\Delta x$ says that the numerical domain of dependence must contain the physical one. It is a much milder restriction than the $\Delta t\propto\Delta x^2$ of the explicit heat equation, and it is enforced ruthlessly: two percent above the limit, an instability seeded by round-off alone overwhelms the solution within a couple of hundred steps.
Its accuracy improves toward the stability limit, and is exact at it. At $C=1$ the update reduces to d'Alembert's formula on the grid and the numerical solution is the exact one. This is a property of the constant-coefficient one-dimensional problem and does not survive variable wave speed, several space dimensions, or irregular meshes, where $C=1$ can no longer hold everywhere at once. But the tendency does survive: for wave propagation, the accurate place to run an explicit scheme is close to the stability limit, not far below it.
Taken together these say something about how a numerical scheme should be judged. The scheme is not accurate because it has a small step, and it is not stable because it is cautious. It is accurate and stable exactly when its discrete structure matches the structure of the equation it is solving, which for a wave equation means moving information along the characteristics at the right speed.
What each lab was about
| Lab | Problem | Method | The point it makes |
|---|---|---|---|
| 1. Nonlinear root finding | Operating point of a diode circuit | Bisection, Newton-Raphson, secant | Robustness and speed are a trade-off; the fastest method needs the most information and the best starting guess |
| 2. Numerical integration | Mean-square and RMS value of a signal | Trapezoidal and Simpson rules | A higher-order rule is not automatically the more accurate one; the structure of the integrand decides |
| 3. ODE numerical solution | Forced vibration of a 1-DOF oscillator | Explicit Euler | A stability limit divides step sizes into those for which the solution converges and those for which it diverges |
| 4. 1D wave equation | Vibration of a stretched string | Central differences in time and space | Accuracy can improve as the stability limit is approached, so a smaller step is not always better |
Key takeaways
- A numerical method is chosen against the information available, not by its convergence rate alone. Bisection needs only a sign change and cannot fail; Newton-Raphson needs a derivative and a good starting guess and repays them with rapid convergence; the secant method buys most of that speed with no derivative at all. All three located the diode operating point at $V_D\approx0.637$ V.
- A higher order does not guarantee a smaller error. On the periodic integrand of Lab 2 the trapezoidal rule was exact as soon as the grid resolved the highest harmonic present, while Simpson's rule — a combination of two trapezoidal results on grids of different spacing — was limited by the coarser of the two and was wrong by more than one unit on the grid where the trapezoidal rule was exact.
- Stability is decided by the eigenvalues and the step size, and nothing else. The explicit Euler oscillator of Lab 3 converged at $\Delta t=0.01$ s and diverged at $\Delta t=0.12$ s, on either side of the limit $\Delta t<0.1$ s that follows from $\lvert1+\lambda\Delta t\rvert<1$. Below the limit the error shrinks in proportion to the step; above it the computed transient grows where the true one decays.
- A stability limit is a threshold, not a guideline. The CFL condition $C\le1$ of Lab 4 is enforced ruthlessly: two percent above it, an instability seeded by round-off alone overwhelms the solution within a couple of hundred steps.
- An exact reference is worth the effort of deriving it. Every lab measured its error against an analytical solution, so that every discrepancy reported could be attributed to the numerical method and to nothing else. Without the exact value $M_x=9$, the wrong Simpson result of Lab 2 would have looked entirely plausible.
- The quality of an answer is decided by how well the discrete method matches the problem, not by how small the step is. This is the one point no single lab makes alone. It is why a periodic integrand suits the trapezoidal rule, why a lightly damped oscillator is awkward for explicit Euler, and why the wave scheme is exact at $C=1$, where its update reduces to d'Alembert's formula on the grid.
Looking ahead. Every method in these labs was written out by hand, and every one of them needed a step size, a stability limit or a convergence check chosen for the problem at hand. The notebooks that follow turn to neural networks, which approximate a solution by a different route: Embedding Physical Knowledge into Neural Networks shows how the physical constraints met here, positivity, conservation, boundary values, can be built into a network, and the physics-informed networks after it put the differential equations of Labs 3 and 4 directly into the training loss, with the finite difference and time-stepping schemes of this notebook serving as the reference they are judged against.