Introduction to Numerical Methods

The vocabulary of this article is the daily currency of computational physics. In molecular dynamics and celestial mechanics, symplectic integrators such as velocity-Verlet and leapfrog beat generic Runge–Kutta schemes because they bound the energy error over billions of steps — a stability property, not an accuracy one, that lets a protein or a planetary system run for astronomical times without drift. In lattice QCD and many-body physics, bound-state energies come from a matrix eigenvalue problem whose conditioning sets how many digits survive, benchmarked against WKB and Bohr–Sommerfeld asymptotics. High-dimensional integrals fall to Monte Carlo, whose $N^{-1/2}$ error is dimension-independent. And every finite-element or finite-difference solver lands on a linear system $Ax=b$ whose condition number $\kappa(A)$ dictates accuracy; the whole industry of preconditioning exists to reshape $A$ toward $\kappa\approx1$ so the iterative solvers behind exascale simulations converge at all. The CFL condition still sets the time step of every explicit weather and astrophysics code on Earth.

Think of a numerical method as photographing a smooth curve onto a finite grid of pixels. Coarse pixels — a large step $h$ — blur the curve, and that blur is truncation error. But every pixel value is also recorded to limited precision, and that is round-off. Crank the resolution up without limit and the per-pixel noise, divided by an ever-smaller $h$, comes to dominate. There is a sweet spot, and the first figure is a portrait of it: two error branches, one falling and one rising, crossing at the optimal step where a one-sided derivative reaches its best possible accuracy of roughly $\sqrt{\varepsilon_{\mathrm{mach}}}\approx10^{-8}$ — half the sixteen digits of double precision, gone before you start.

Stability is a different question entirely — it is about memory. In a time-stepping calculation, does a mistake made early get forgotten or amplified? A stable scheme behaves like a damped system that forgets perturbations; an unstable one like an inverted pendulum where the tiniest nudge explodes. The cleanest demonstration is forward Euler on a decaying exponential: below a critical step size the numerical solution tracks the true decay, but one step too large and it flips into a growing oscillation — a divergence with nothing to do with formal accuracy and everything to do with the amplification factor per step. The second figure shows both sides of that knife-edge, with the exact solution for reference.

Conditioning, finally, lives in the problem and not the solver. Subtracting $1.0000001$ from $1.0000002$ loses six of seven significant figures before a line of code runs — catastrophic cancellation, which no algorithm recovers from; the only cure is to reformulate the problem so it never happens.

Numerical differentiation is the cleanest place to watch truncation error appear. The simplest estimate of a derivative replaces the limit in its definition with a small finite step, and Taylor's theorem tells us exactly what we lose.

where: $f'(x)$ — the true derivative; $h$ — the step size, a small positive number; $O(h)$ — the leading error, proportional to $\tfrac12 h\,f''$. Reading it: the forward difference is first order, so halving $h$ only halves the error — a poor return a symmetric stencil improves for free. Sampling the function on both sides of $x$ cancels the leading error term by symmetry, buying an extra order at the same two function evaluations. This is the workhorse first-derivative formula.

where: $O(h^2)$ — leading error, proportional to $\tfrac16 h^2 f'''$. Reading it: the central difference is second order, so halving $h$ quarters the error — one order higher than (1) at no extra cost. Every discretization carries a truncation term that falls with $h$ and, in finite precision, a round-off term that grows as $h$ shrinks. Writing both explicitly gives the most important relation in the subject — the error budget.

where: $C,D$ — problem-dependent constants of order unity; $p$ — the order of accuracy; $\varepsilon_{\mathrm{mach}}$ — machine epsilon, the relative round-off unit. Reading it: the first term is the blur, the second the noise; their sum has a genuine minimum at finite $h$, so "smaller $h$" is not a strategy but a trap. The round-off unit is a hard property of the number format, not a fudge factor. In IEEE 754 double precision it is the gap between $1$ and the next representable number.

where: $2^{-52}$ — the spacing of the 52-bit mantissa at exponent zero. Reading it: this fixes the noise floor of every calculation at about sixteen decimal digits, and through (3) the best accuracy any single differencing operation can reach. Integration discretizes the opposite way — a continuous area becomes a weighted sum of samples. Joining the samples by straight lines gives the composite trapezoidal rule, whose error is set by the curvature it ignores.

where: $h=(b-a)/n$ — uniform panel width; $f_i=f(a+ih)$ — sampled values. Reading it: endpoints get weight $1$, interior points weight $2$; the global error is $-\tfrac{(b-a)h^2}{12}f''(\eta)$, so a straight-line fit is second order. Fitting parabolas through triples of points cancels two more orders, giving Simpson's rule — the best return-on-effort quadrature for smooth integrands.

where: $n$ — the number of panels, which must be even; the $1,4,2,4,\dots,4,1$ weights are the parabola quadrature coefficients. Reading it: the error is $O(h^4)$, so halving $h$ divides it by sixteen — two orders better than the trapezoid for one extra line of code. Differential equations are marched forward in time by turning the derivative into a difference and stepping. The explicit (forward) Euler method is the prototype of every such integrator.

where: $y_n\approx y(t_n)$ — the numerical solution at time $t_n$; $f$ — the right-hand side of $y'=f(t,y)$. Reading it: the local per-step error is $O(h^2)$, but it accumulates over $N=(T-t_0)/h$ steps into a global error $O(h)$ — Euler is first order, and its stability is a separate matter. Whether the problem amplifies error is a property of the linear system, not the algorithm. For $Ax=b$ that number is the condition number.

where: $\lVert\cdot\rVert$ — an induced matrix norm, often the 2-norm; $\kappa\ge1$ always, larger meaning more ill-conditioned. Reading it: $\kappa$ is the worst-case amplification of relative input error, so $\kappa=10^{8}$ can cost eight digits however carefully the system is solved. For time-dependent partial differential equations, stability couples the time step to the spatial grid. The Courant–Friedrichs–Lewy condition is the archetypal constraint for an explicit wave or advection scheme.

where: $c$ — the wave or advection speed; $\Delta t,\Delta x$ — the time and space step sizes. Reading it: information must not physically outrun the stencil in one step, so refining the spatial grid forces a proportionally smaller time step.

The behaviour of the total error as the step size or precision changes lays bare the central trade-offs, and the condition number sets the ultimate limit that no algorithm escapes.

1 · Truncation, order, and the round-off floor (confidence builder)

Estimate $f'(1)$ for $f(x)=e^{x}$ (exact value $e=2.718282$) by finite differences, and show the forward and central formulas obey their advertised orders.

Step 1. Forward difference at $h=0.1$: $\dfrac{e^{1.1}-e^{1.0}}{0.1}=\dfrac{3.004166-2.718282}{0.1}=2.858842$, an error of $+0.140560$. The predicted truncation error is $\tfrac12 h f''=\tfrac12(0.1)(2.718282)=0.135914$ — the same size, confirming the $O(h)$ term dominates. Forward diff error ≈ 0.14.

Step 2. Central difference at $h=0.1$: $\dfrac{e^{1.1}-e^{0.9}}{0.2}=\dfrac{3.004166-2.459603}{0.2}=2.722815$, an error of $+0.004533$. The predicted $\tfrac16 h^2 f'''=\tfrac16(0.01)(2.718282)=0.004530$ matches to three figures — and the error is 31 times smaller than the forward difference at the same two evaluations. Central diff error ≈ 0.0045.

Step 3. Halve nothing — shrink $h$ by ten. Central difference at $h=0.01$: $\dfrac{e^{1.01}-e^{0.99}}{0.02}=2.718327$, an error of $+4.53\times10^{-5}$. The error fell by a factor of almost exactly $100=10^{2}$ when $h$ fell by $10$ — the signature of second order. Cutting $h$ tenfold cut the central-difference error 100-fold, as $O(h^2)$ demands.

Step 4. Do not chase this forever. By equation (3) the round-off branch $\varepsilon_{\mathrm{mach}}/h$ turns the total error around near $h_{\mathrm{opt}}\sim\sqrt{\varepsilon_{\mathrm{mach}}}\approx1.5\times10^{-8}$ for the first derivative; below that step the difference of two near-equal exponentials is pure noise and the estimate gets worse. The right way to more digits is a higher-order stencil, not a smaller step.

2 · Trapezoid versus Simpson on a known integral

Evaluate $\displaystyle\int_0^1 e^{x}\,\mathrm{d}x$ (exact value $e-1=1.718282$) with the trapezoidal and Simpson rules, and verify their orders by refining the grid.

Step 1. Trapezoid with $n=2$ panels ($h=0.5$): $\dfrac{h}{2}\big[f_0+2f_1+f_2\big]=\dfrac{0.5}{2}\big[1+2(1.648721)+2.718282\big]=1.753931$, an error of $+0.035649$. Trapezoid, coarse: error ≈ 3.6×10⁻².

Step 2. Refine to $n=4$ ($h=0.25$): $\dfrac{0.25}{2}\big[f_0+2(f_1+f_2+f_3)+f_4\big]=1.727222$, an error of $+0.008940$. The ratio of errors is $0.035649/0.008940=3.99\approx4=2^{2}$ — halving $h$ quartered the error, exactly the $O(h^2)$ law. Trapezoid halving ÷ 4: confirmed second order.

Step 3. Simpson with $n=2$ ($h=0.5$): $\dfrac{h}{3}\big[f_0+4f_1+f_2\big]=\dfrac{0.5}{3}\big[1+4(1.648721)+2.718282\big]=1.718861$, an error of $+5.79\times10^{-4}$ — already 60 times better than the trapezoid on the same three points. Simpson, coarse: error ≈ 5.8×10⁻⁴.

Step 4. Refine Simpson to $n=4$ ($h=0.25$): $1.718319$, an error of $+3.70\times10^{-5}$. The ratio $5.79\times10^{-4}/3.70\times10^{-5}=15.7\approx16=2^{4}$ — halving $h$ cut the error sixteen-fold, the $O(h^4)$ law in action. Simpson halving ÷ 16: confirmed fourth order — two free orders for one $1,4,2,4,1$ weighting.

Numerical analysis had many founders — Newton and Euler gave the first iterations and integrators, Gauss the first quadratures and elimination — but the theory of error this article is really about had one: the man who first explained, precisely, why some algorithms return trustworthy answers to hard problems and others do not.

Texts. Press, Teukolsky, Vetterling & Flannery, Numerical Recipes, 3rd ed. (Chs. 5, 17, 20) · Trefethen & Bau, Numerical Linear Algebra (SIAM 1997) · Stone & Goldbart, Mathematics for Physics (CUP 2009) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed. · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed. (Ch. 27). Key numerical benchmarks: machine epsilon $2^{-52}=2.220\times10^{-16}$ and $\sqrt{\varepsilon_{\mathrm{mach}}}=1.490\times10^{-8}$; central-difference error $O(h^2)$ and trapezoid/Simpson errors $O(h^2)/O(h^4)$ from Taylor remainders; worked-example values $2.858842$, $2.722815$, $2.718327$ for $f'(1)$ of $e^x$, and $1.753931$, $1.727222$, $1.718861$, $1.718319$ for $\int_0^1 e^x\,\mathrm{d}x$ (exact $1.718282$), with error ratios $3.99$ and $15.7$; the forward-Euler stability bound $h\le2/|\lambda|$ from $|1+h\lambda|\le1$; and Vieta's relation $x_1x_2=c/a$ for the stable quadratic formula.

Prerequisites. Order Symbols and Gauge Functions · Taylor's theorem · Linear algebra. Leads to. Asymptotic Series and Poincaré Expansions · Optimal Truncation and Superasymptotics · WKB Approximation.

Ask the tutor: ask me to walk through the error-budget curve for any specific derivative or integral, and we will find its optimal step size and best achievable accuracy together.