Numerical Approximation and Error

Sources of numerical error: truncation, floating-point roundoff, stability, and convergence in computation.

What it is

Numerical approximation is the physical and mathematical foundation that models, bounds, and manages the unavoidable discrepancies arising when continuous physical laws formulated over $\mathbb{R}$ are evaluated on finite-precision computing architectures.

The total computational error splits into three irreducible sources:

1. Truncation error ($E_{\text{trunc}}$): Arises from replacing continuous infinite operators (infinite Taylor series, derivatives as limits $\lim_{h \to 0}$, or Riemann integrals) with finite discrete approximations parameterized by mesh step $h$. Its asymptotic decay follows Landau notation $\mathcal{O}(h^p)$, where $p$ is the convergence order.
2. Roundoff error ($E_{\text{round}}$): Stemming from representing real continuum numbers using fixed binary formats under IEEE-754 standard (single precision Float32 with 24-bit significand and machine epsilon $\epsilon_{\text{mach}} \approx 1.19 \times 10^{-7}$; or double precision Float64 with 53-bit significand and $\epsilon_{\text{mach}} \approx 2.22 \times 10^{-16}$). Every arithmetic operation incurs relative roundoff bounded by $\epsilon_{\text{mach}}$.
3. Catastrophic cancellation: A pathological loss of precision occurring when subtracting two nearly equal floating-point numbers: identical leading significant digits cancel, leaving only prior roundoff noise in the remaining bits.

Why it appears in physics

In modern computational physics (astrophysics, quantum many-body systems, climate modeling, fusion plasma physics), nonlinear differential systems lack closed-form analytical solutions:

• Simulating $N$-body celestial orbits, turbulence, or tokamak magnetic confinement demands spatial and temporal discretization.
• A fundamental tension exists known as the step-size trade-off: reducing step size $h$ drives truncation error down as $C_1 h^p$, but arithmetic operation count scales as $1/h$, compounding floating-point roundoff error as $C_2 \epsilon_{\text{mach}} / h^q$. Below an optimal step $h^*$, shrinking $h$ exacerbates total numerical error.
• Neglecting numerical stability and matrix condition numbers leads to unphysical simulations characterized by energy divergence, violated conservation laws, or artificial numerical blowup.

Visual or geometric idea

In a log-log plot ($\log E$ vs $\log h$), total error forms a characteristic "V"-shaped profile:

• At large $h$, truncation error dominates along a line with positive slope $+p$: $\log E \approx p \log h + \text{const}$.
• At tiny $h$, floating-point cancellation and roundoff compound, creating a branch with negative slope.
• The global minimum marks the optimal step size $h^*$. Stepping below $h^*$ consumes CPU cycles while degrading numerical accuracy.

Minimal symbolic language

1. Absolute and relative error: $E_{\text{abs}} = |x_{\text{true}} - x_{\text{approx}}|, \qquad E_{\text{rel}} = \frac{|x_{\text{true}} - x_{\text{approx}}|}{|x_{\text{true}}|}$
2. IEEE-754 standard floating-point model: $\operatorname{fl}(x) = x (1 + \delta), \qquad |\delta| \leq \epsilon_{\text{mach}}$
3. Catastrophic cancellation ($x \approx y$): $\operatorname{fl}(x - y) = (x - y) + x \delta_1 - y \delta_2 = (x - y) \left[ 1 + \frac{x \delta_1 - y \delta_2}{x - y} \right]$ As $(x - y) \to 0$, relative error diverges unless expressions are reformulated algebraically.
4. Condition number of problem $f(x)$: $\kappa(x) = \left| \frac{x f'(x)}{f(x)} \right|$ If $\kappa(x) \gg 1$, the physical problem is ill-conditioned (hypersensitive to perturbations).

What the result means

Understanding numerical error elevates scientific simulations from ad-hoc coding into verifiable science:

• A computational calculation possesses no scientific validity without an accompanying convergence proof or quantified error bound.
• In orbital or molecular dynamics, confirming that global error drops at the theoretical rate ($h^2, h^4$) validates both software implementation and numerical stability.

When not to use it

1. Do not blindly decrease $h$ past the roundoff threshold: In numerical finite-difference differentiation, selecting $h = 10^{-16}$ in Float64 yields zero or random noise due to catastrophic cancellation.
2. Do not use single precision (Float32) for long-duration dynamic trajectories: Roundoff drift accumulated over millions of integration steps destroys energy conservation in physical orbits.
3. Do not conflate numerical error with physics model error: Numerical error quantifies the gap between algorithm and differential equation; it cannot detect missing friction or neglected quantum corrections in the original physics formulation.

Operational summary

• Write the continuous differential or integral equation and choose an approximation scheme.
• Identify the asymptotic truncation order $\mathcal{O}(h^p)$.
• Inspect formulas for catastrophic cancellation hazards (differences of close numbers, square root subtractions); eliminate them via algebraic conjugates or local Taylor series.
• Choose step size $h$ balancing truncation against IEEE-754 machine epsilon limits.
• Perform a systematic grid refinement test ($h, h/2, h/4$) to experimentally confirm expected convergence orders.

Procedure for Numerical Error Control and Mitigation

Before calculating

Before executing computational physical simulations:

• Identify runtime floating-point precision (Float32: $\epsilon_{\text{mach}} \approx 1.2 \times 10^{-7}$, 7 decimal digits; Float64: $\epsilon_{\text{mach}} \approx 2.2 \times 10^{-16}$, 16 decimal digits).
• Scan mathematical formulas for catastrophic cancellation vulnerabilities (e.g. $1 - \cos(x)$ for $x \to 0$, or $\sqrt{x + 1} - \sqrt{x}$ for $x \gg 1$).

Step 1: identify the mathematical object

Classify the numerical algorithm and its truncation profile:

• Finite-difference derivative: Truncation $\mathcal{O}(h)$ for forward/backward or $\mathcal{O}(h^2)$ for centered differences.
• Numerical quadrature: Truncation $\mathcal{O}(h^2)$ for trapezoidal or $\mathcal{O}(h^4)$ for Simpson's rule.
• ODE solver (Runge-Kutta): Local truncation error $\mathcal{O}(h^{p+1})$, global error $\mathcal{O}(h^p)$.

Step 2: choose the operation (Cancellation elimination)

If nearly equal values are subtracted, algebraically reformulate before coding:

• Root subtraction: $\sqrt{a} - \sqrt{b} = \frac{a - b}{\sqrt{a} + \sqrt{b}}$ (multiplying by conjugate removes cancellation).
• Small-angle trigonometry: $1 - \cos(x) = 2 \sin^2(x/2)$ or Taylor expansion $\frac{x^2}{2} - \frac{x^4}{24} + \dots$.
• Quadratic formula with $b^2 \gg 4 a c$: $x_1 = \frac{-b - \operatorname{sgn}(b)\sqrt{b^2 - 4ac}}{2a}, \qquad x_2 = \frac{c}{a x_1}$

Step 3: preserve units and signs

• Absolute error $E_{\text{abs}}$ preserves identical dimensions of the physical observable.
• Relative error $E_{\text{rel}}$ and machine epsilon $\epsilon_{\text{mach}}$ are dimensionless pure numbers.
• In numerical differentiation $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$, step size $h$ must carry identical physical units as input variable $x$.

Step 4: calculate (Optimal step-size calculation $h^*$)

For a numerical method with truncation error $\alpha h^p$ and roundoff growth $\frac{\epsilon_{\text{mach}}}{h^q}$:

1. Formulate total error: $E_{\text{total}}(h) = \alpha h^p + \beta \frac{\epsilon_{\text{mach}}}{h^q}$
2. Set derivative to zero: $\frac{d E_{\text{total}}}{dh} = p \alpha h^{p-1} - q \beta \frac{\epsilon_{\text{mach}}}{h^{q+1}} = 0$
3. Solve for optimal step size: $h^* = \left( \frac{q \beta \epsilon_{\text{mach}}}{p \alpha} \right)^{\frac{1}{p + q}}$ (For centered differences with $p=2, q=1$: $h^* \sim \epsilon_{\text{mach}}^{1/3}$. In Float64, $h^* \sim (10^{-16})^{1/3} \approx 10^{-5} - 10^{-6}$).

Step 5: interpret the result (Richardson convergence check)

Compute approximations across three halved step sizes ($h, h/2, h/4$) denoted $A_1, A_2, A_3$:

1. Calculate convergence ratio: $R = \frac{A_1 - A_2}{A_2 - A_3} \approx 2^p$
2. If $R \approx 2^p$ (e.g. $R \approx 4$ for 2nd-order $p=2$, or $R \approx 16$ for 4th-order $p=4$), the computation resides cleanly in the asymptotic convergence regime free from roundoff contamination.
3. Richardson extrapolation for leading error cancellation: $A_{\text{refined}} = A_3 + \frac{A_3 - A_2}{2^p - 1}$

Self-check checklist

• Is selected step $h$ larger than the roundoff minimum $h^*$?
• Did I eliminate catastrophic cancellations via conjugate or series expansions?
• Does empirical convergence ratio $R$ match theoretical $2^p$?
• Did I enforce Float64 precision for long dynamic runs?
• Are physical conservation invariants (energy, momentum) preserved within discretization bounds?

Guided Example: Catastrophic Cancellation in Relativistic Kinematics

Physical Situation

A subatomic particle of rest mass $m = 1.000 \text{ GeV}/c^2$ is accelerated inside the Large Hadron Collider (LHC) to total relativistic energy $E = 1000.000000000100 \text{ GeV}$. We compute relativistic momentum $p$ using Einstein's energy-momentum relation: $E^2 = p^2 c^2 + m^2 c^4 \implies p c = \sqrt{E^2 - m^2 c^4}$ In natural units ($c=1$), $p = \sqrt{E^2 - m^2}$. Because $E \gg m$, momentum $p$ is nearly indistinguishable from $E$. If kinetic energy and velocity deficit $\Delta v = c - v = 1 - \frac{p}{E} = 1 - \sqrt{1 - \frac{m^2}{E^2}}$ are computed naively in floating-point software, catastrophic cancellation destroys precision.

Mathematical Translation

Evaluating with $E = 10^4 \text{ GeV}$ and $m = 1 \text{ GeV}$: $\gamma = \frac{E}{m} = 10^4, \qquad \frac{m^2}{E^2} = 10^{-8}$ The naive subtraction evaluates: $1 - \sqrt{1 - 10^{-8}}$ In floating-point math, $\sqrt{1 - 10^{-8}} = 0.9999999950000000\dots$ Subtracting $1 - 0.999999995\dots$ cancels the leading 8 identical digits, destroying half the mantissa in Float64 and producing total cancellation to $0.0$ in Float32.

Step-by-Step Calculation

1. Conjugate Algebraic Reformulation: $\Delta v = 1 - \sqrt{1 - \frac{m^2}{E^2}} = \frac{\left(1 - \sqrt{1 - \frac{m^2}{E^2}}\right)\left(1 + \sqrt{1 - \frac{m^2}{E^2}}\right)}{1 + \sqrt{1 - \frac{m^2}{E^2}}} = \frac{m^2 / E^2}{1 + \sqrt{1 - \frac{m^2}{E^2}}}$
2. Stable Floating-Point Evaluation: The denominator is now a stable sum of positive quantities: $1 + \sqrt{1 - 10^{-8}} \approx 1 + 0.999999995 = 1.999999995$ The numerator is: $\frac{m^2}{E^2} = 10^{-8}$ Therefore: $\Delta v = \frac{10^{-8}}{1.999999995} \approx 5.0000000125 \times 10^{-9} \text{ (in units of } c)$
3. Binomial Taylor Alternative: $1 - \sqrt{1 - \epsilon} = \frac{1}{2}\epsilon + \frac{1}{8}\epsilon^2 + \dots = \frac{1}{2}(10^{-8}) + \frac{1}{8}(10^{-16}) = 5.0000000000000125 \times 10^{-9}$

Reading the Result

The conjugate reformulation preserves all 16 significant decimal figures in Float64. The particle travels just $\Delta v \approx 5.0 \times 10^{-9} c \approx 1.5 \text{ m/s}$ below the speed of light. Naive evaluation would have returned zero or suffered $50\%$ rounding error.

Common Error Avoided

Directly evaluating differences of square roots of nearly identical quantities. Rationalizing with algebraic conjugates eradicates catastrophic cancellation at zero additional computational cost.

Inverse Example: Identifying the Optimal Step Size in Numerical Differentiation

Physical Situation

A spacecraft flight computer estimates acceleration $a(t) = v'(t)$ from velocity telemetry $v(t) = \sin(t)$ at $t = 1.0 \text{ s}$ using central differences on a standard 64-bit IEEE 754 platform ($\epsilon_{\text{mach}} \approx 2.22 \times 10^{-16}$). The exact analytical derivative is $v'(1.0) = \cos(1.0) \approx 0.5403023058681398$.

Given Result

Central difference formula: $v'(t) \approx \frac{v(t + h) - v(t - h)}{2h}$ Testing three step sizes:

• $h = 10^{-2}$: Absolute error $\epsilon = 9.0 \times 10^{-6}$ (dominated by $\mathcal{O}(h^2)$ truncation).
• $h = 10^{-5}$: Absolute error $\epsilon = 4.2 \times 10^{-12}$ (optimal balance).
• $h = 10^{-14}$: Absolute error $\epsilon = 0.08$ (dominated by $\epsilon_{\text{mach}} / h$ floating-point cancellation).

Analytical Reconstruction

Total numerical error balances truncation and roundoff: $E(h) = \frac{M h^2}{6} + \frac{\epsilon_{\text{mach}}}{h}$ Minimizing $E'(h) = \frac{M h}{3} - \frac{\epsilon_{\text{mach}}}{h^2} = 0$: $h_{\text{opt}} = \left( \frac{3 \epsilon_{\text{mach}}}{M} \right)^{1/3} \sim \epsilon_{\text{mach}}^{1/3} \approx (10^{-16})^{1/3} \approx 10^{-5.3}$

What Physical Question It Can Answer

Why does making the numerical step arbitrarily small ($h \to 0$) destroy the calculated derivative? Because dividing rounding noise by an infinitesimal $h$ amplifies machine noise into macroscopic errors, as proved by the $h = 10^{-14}$ calculation.

Boundary Verification

As $h \to 0$, $E(h) \to \infty$; as $h \to 1$, $E(h) \sim h^2 \to \mathcal{O}(1)$. The minimum lies strictly at $h_{\text{opt}} \sim \epsilon_{\text{mach}}^{1/3}$.

Common Mistakes in Numerical Approximation and Error

1. Shrinking step size $h$ assuming error monotonically vanishes

It seems reasonable because elementary calculus defines derivatives and integrals as the limit $h \to 0$. It fails because on real finite-precision machines (IEEE-754), total computational error traces a "V"-shape: as $h$ decreases, truncation shrinks but roundoff error compounds as $\epsilon_{\text{mach}} / h^q$. Below optimal step $h^*$, smaller steps amplify global error. It is detected when empirical error stops decreasing and begins oscillating or rising steeply. It is corrected by determining optimal step size $h^* \sim \epsilon_{\text{mach}}^{1/(p+q)}$ and running grid convergence checks. Mini-example: evaluating $f'(x) \approx \frac{f(x+h) - f(x)}{h}$ with $h = 10^{-16}$ in Float64 yields $0.0$ instead of the true slope, because $f(x+10^{-16})$ rounds back to $f(x)$.

2. Subtracting nearly equal floating-point numbers (catastrophic cancellation)

It seems reasonable because computers are deterministic machines providing 16 decimal digits of precision. It fails because if two floats match in their first 14 decimal digits, subtracting them cancels those 14 digits, leaving only 2 digits of genuine precision followed by roundoff noise. It is detected when sudden accuracy loss occurs in root differences, small-angle trigonometry, or quadratic formulas. It is corrected by algebraically reformulating expressions via conjugates, Taylor expansions, or half-angle identities. Mini-example: in Float32, evaluating $1 - \cos(10^{-4})$ outputs $0.0$; computing $2 \sin^2(10^{-4}/2)$ yields accurate $5.000000 \times 10^{-9}$.

3. Testing floating-point equality with exact comparisons

It seems reasonable because mathematically equivalent expressions should evaluate identically. It fails because binary floating-point arithmetic is neither associative nor distributive: $(0.1 + 0.2) + 0.3 \neq 0.1 + (0.2 + 0.3)$ under IEEE-754. Evaluating 0.1 + 0.2 == 0.3 returns false in most programming languages. It is detected when iterative convergence loops hang in infinite cycles. It is corrected by using absolute and relative error bounds: abs(x - y) <= atol + rtol * abs(y). Mini-example: an iteration while x != 1.0 with x += 0.1 never terminates because x never hits binary 1.0000000000000000.

4. Confusing numerical stability with approximation order

It seems reasonable because higher-order solvers (like RK4) are more accurate than low-order schemes (like Euler). It fails because an algorithm can possess high asymptotic accuracy while being unstable for stiff problems or if step size exceeds the algorithm's stability boundary. It is detected when numerical traces explode abruptly to $\pm \infty$ after behaving normally in initial steps. It is corrected by inspecting the linear stability domain $z = \lambda \Delta t$ in the complex plane and adopting implicit schemes (Crank-Nicolson, Radau) for stiff equations. Mini-example: in chemical kinetics with rate constants $k_1 = 1, k_2 = 10^6$, explicit RK4 requires $\Delta t < 10^{-6}$ not for accuracy, but to prevent explosive numerical blowup.

5. Summing series in descending order of magnitude

It seems reasonable because adding dominant large terms first feels like capturing the bulk of the result. It fails because adding a tiny term to a massive accumulator shifts the small term's bits off the significand, causing complete numerical absorption. It is detected when summing $N = 10^7$ small positive floats causes the accumulator to stall and stop growing. It is corrected by sorting terms in ascending order (smallest to largest) or implementing Kahan compensated summation. Mini-example: in Float32, adding $1.0$ to an accumulator of $10^8$ has zero effect: $10^8 + 1.0 = 10^8$. Adding $1.0$ one million times in a loop still leaves the accumulator at $10^8$.

Quick detection rule

• Subtracting two values? Check if $|x/y| \in [0.99, 1.01]$. If so, catastrophic cancellation is occurring.
• Reducing step size $h$? Plot error on a log-log axis. If the line turns upward, you have breached the roundoff floor.

Physical Applications of Numerical Approximation and Error

1. Astrodynamics: interplanetary trajectory calculations and orbital semimajor axes

Appears in aerospace engineering when integrating equations of motion for deep-space probes (Voyager, New Horizons) under the gravitational influence of the Sun and Solar System planets. It is calculated by integrating the $N$-body differential equations $\ddot{\mathbf{r}} = -\sum G M_i \frac{\mathbf{r} - \mathbf{r}_i}{|\mathbf{r} - \mathbf{r}_i|^3}$ over multi-decade mission horizons. The result provides probe state vectors for gravity-assist flybys with trajectory error bounds within tens of kilometers after traversing billions of kilometers. The critical danger is accumulated roundoff and artificial energy dissipation caused by non-conservative integrators (standard RK4); over a 30-year mission with hourly steps, tiny $10^{-15}$ roundoff per step compounds into noticeable orbital energy drift. As a final sanity check, high-order symplectic integrators (Gauss-Radau, Yoshida) are deployed and strict invariance of the Laplace-Runge-Lenz vector and total Hamiltonian energy is enforced.

2. Computational quantum mechanics: stationary Schrödinger eigenvalue problems

Appears in atomic and solid-state physics when solving for energy eigenstates and wavefunctions $\psi(x)$ in non-harmonic potentials or diatomic molecules via the shooting method or spatial finite differences: $-\frac{\hbar^2}{2m} \frac{d^2\psi}{dx^2} + V(x)\psi = E\psi$. It is calculated using the sixth-order Numerov algorithm $\mathcal{O}(h^6)$ for second-order differential equations lacking first derivatives. The result determines discrete bound-state energies $E_n$ and quantum barrier tunneling probabilities in semiconductor heterostructures. The critical danger is using too coarse a spatial grid $h$ in classically forbidden regions where wavefunctions decay exponentially; if truncation error fails to resolve exponential curvature, the numerical solution blows up toward $\pm \infty$, preventing eigenvalue convergence. As a final sanity check, wavefunction orthonormality $\int |\psi_n|^2 dx = 1$ is verified and numerical eigenvalues are benchmarked against first-order perturbation theory.

3. Geophysical fluid dynamics: numerical weather prediction and climate models

Appears in atmospheric physics when solving continuous hydrodynamic Navier-Stokes primitive equations across three-dimensional global atmospheric grids. It is calculated by discretizing pressure gradient and wind advection terms via centered finite differences or spherical harmonic spectral transforms. The result predicts the temporal evolution of storm tracks, thermal fronts, and cyclonic systems across multi-day forecast windows. The critical danger is violating the Courant-Friedrichs-Lewy condition (CFL: $C = \frac{u \Delta t}{\Delta x} \leq 1$); if wind or acoustic/gravity wave speeds exceed one grid cell per time step, truncation errors amplify exponentially, triggering explosive unphysical numerical blowup. As a final sanity check, adaptive time-stepping tracks local CFL numbers and hyperdiffusion operators selectively damp sub-grid eddies without altering planetary Rossby wave dynamics.

4. Plasma physics and fusion: charged particle orbit tracking in Tokamak reactors

Appears in magnetic confinement fusion (ITER) when calculating guiding-center and full-orbit trajectories of thermal deuterons and fast fusion alpha particles trapped in helical toroidal fields. It is calculated by integrating the Lorentz force $\mathbf{F} = q (\mathbf{E} + \mathbf{v} \times \mathbf{B})$ via the Boris kinematic algorithm. The result resolves gyromotion Larmor radii and predicts whether trapped banana orbits drift into divertor or first-wall armor tiles (confinement loss). The critical danger is employing standard Runge-Kutta solvers for the rotational term $\mathbf{v} \times \mathbf{B}$; standard solvers fail to conserve magnetic moment $\mu = m v_\perp^2 / (2B)$ and gyro-radius, causing particles to gain artificial numerical kinetic energy and escape confinement. As a final sanity check, Boris integration exactly preserves velocity magnitude in pure magnetic fields ($|\mathbf{v}_{n+1}| = |\mathbf{v}_n|$), guaranteeing strict kinetic energy conservation across millions of gyro-periods.

5. Electromagnetic optics: FDTD modeling of photonic metamaterials and waveguides

Appears in photonics and optical communications when simulating laser pulse propagation across dielectric waveguides and photonic bandgap crystals via Yee-grid Finite-Difference Time-Domain (FDTD). It is calculated by interleaving spatial and temporal components of electric $\mathbf{E}$ and magnetic $\mathbf{H}$ fields on staggered cubic grids shifted by half steps in space and time. The result provides spectral transmission curves, dispersion relations, and chromatic group velocity dispersion in specialty optical fibers. The critical danger is numerical dispersion: on discrete grids, the phase velocity of simulated electromagnetic waves $v_{\text{num}}$ does not equal $c$ identically, but varies artificially with propagation angle and wavelength-to-grid ratio $\lambda / \Delta x$. On coarse grids, ultra-short pulses artificially broaden and distort due to spatial discretization truncation. As a final sanity check, spatial grid resolution is constrained to at least 20 cells per dielectric wavelength ($\Delta x \leq \lambda_{\text{min}} / (20 n_{\text{ref}})$) and Courant numbers are locked to the multidimensional stability optimum.

Guided Practice: Numerical Approximation and Error

Level 1: recognition

Classify the dominant numerical error in each computational operation (truncation, roundoff, or catastrophic cancellation):

1. Computing $\pi$ by summing the first 100 terms of the Gregory-Leibniz alternating series $\sum_{k=0}^{99} \frac{(-1)^k}{2k+1}$.
2. Evaluating $y = \sqrt{1000001} - 1000$ on a standard computer.
3. Storing $f(x) = \frac{1}{3} x$ as Float32 in memory.
4. Integrating a Newtonian planetary orbit over 10 million cycles observing secular energy drift.

Solution:

1. Truncation error (terminating an infinite series after finite operations).
2. Catastrophic cancellation (subtracting two nearly identical large floating-point numbers).
3. Roundoff error (the fraction $1/3$ is an infinite binary recurring float $0.010101\dots_2$ truncated to 24 significand bits).
4. Accumulated roundoff drift (random-walk and systematic accumulation over millions of floating-point operations).

Level 2: direct application

Approximate the first derivative of $f(x) = e^x$ at $x=0$ via forward finite differences: $D_+(h) = \frac{e^h - 1}{h}$ Exact theoretical derivative is $f'(0) = e^0 = 1$. Derive the analytical truncation error via Taylor expansion of $e^h$ and compute numerical error for $h = 0.1$ and $h = 0.01$.

Solution:

1. Taylor series: $e^h = 1 + h + \frac{h^2}{2} + \frac{h^3}{6} + \mathcal{O}(h^4)$
2. Substitute into operator: $D_+(h) = \frac{(1 + h + \frac{h^2}{2} + \dots) - 1}{h} = 1 + \frac{h}{2} + \mathcal{O}(h^2)$
3. Truncation error: $E_{\text{trunc}} = |D_+(h) - 1| = \frac{h}{2} + \mathcal{O}(h^2)$ First-order accuracy $\mathcal{O}(h)$.
4. Numerical evaluation:
   • For $h = 0.1$: $E_{\text{trunc}} \approx 0.1 / 2 = 0.05$ ($5\%$ error; computed value $1.05171$).
   • For $h = 0.01$: $E_{\text{trunc}} \approx 0.01 / 2 = 0.005$ ($0.5\%$ error; computed value $1.005017$). Decimating step size by 10 cuts error by exactly 10 (linear convergence).

Level 3: tool selection

In astrophysical software, you must evaluate gravitational potential energy differences between close orbital points ($r_2 = r_1 + \delta$ with $\delta \ll r_1$): $\Delta U = -G M m \left( \frac{1}{r_2} - \frac{1}{r_1} \right)$ How should this expression be rewritten to guarantee numerical stability free from cancellation?

Solution: Evaluating $\frac{1}{r_2} - \frac{1}{r_1}$ directly causes catastrophic floating-point cancellation. Algebraically combine fractions before computing: $\frac{1}{r_2} - \frac{1}{r_1} = \frac{r_1 - r_2}{r_1 r_2} = \frac{-\delta}{r_1 (r_1 + \delta)}$ Substituting back: $\Delta U = \frac{G M m \delta}{r_1 (r_1 + \delta)}$ This expression uses only products and additions of positive numbers, eliminating subtraction of close numbers and retaining all 16 Float64 digits, even when $\delta = 10^{-12} r_1$.

Level 4: physical interpretation

A computational physicist simulates a simple harmonic oscillator $m \ddot{x} + k x = 0$ via explicit Euler stepping ($x_{n+1} = x_n + v_n \Delta t$; $v_{n+1} = v_n - \omega^2 x_n \Delta t$). Over 1000 cycles, the oscillation amplitude diverges exponentially and mechanical energy grows by a factor of 10. Explain why this behavior is not floating-point roundoff but algorithmic numerical instability.

Solution: Compute the state transition matrix for explicit Euler: $\begin{pmatrix} x_{n+1} \\ v_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & \Delta t \\ -\omega^2 \Delta t & 1 \end{pmatrix} \begin{pmatrix} x_n \\ v_n \end{pmatrix}$ Matrix determinant is: $\det(A) = 1 \cdot 1 - (\Delta t)(-\omega^2 \Delta t) = 1 + \omega^2 \Delta t^2$ The determinant represents the phase-space volume expansion per step. Because $\det(A) = 1 + \omega^2 \Delta t^2 > 1$ for any finite $\Delta t > 0$, eigenvalues have magnitude strictly greater than 1: $|\lambda| = \sqrt{1 + \omega^2 \Delta t^2} > 1$. Physical interpretation: Explicit Euler injects artificial numerical energy at every step. This is an algorithmic instability, not roundoff. Conservative physical systems require symplectic integrators (Symplectic Euler, Verlet) where $\det(A) \equiv 1$.

Level 5: intentional mistake

To solve $a x^2 + b x + c = 0$ with $a = 1$, $b = 10^8$, $c = 1$, a student codes the standard textbook formula: $x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a}, \qquad x_2 = \frac{-b - \sqrt{b^2 - 4ac}}{2a}$ The code outputs $x_1 = 0.0$ and $x_2 = -10^8$. Prove that $x_1 = 0$ is completely wrong and determine the correct root using a stabilized algorithm.

Conceptual blunder: Because $b = 10^8$, $b^2 = 10^{16}$. In Float64 arithmetic, subtracting $4 a c = 4$ from $10^{16}$ falls beyond the 53-bit significand, so $\sqrt{b^2 - 4ac}$ evaluates to $10^8$ identically. Evaluating $x_1 = \frac{-10^8 + 10^8}{2} = 0.0$ destroys the smaller root via catastrophic cancellation. However, Vieta's theorem requires $x_1 \cdot x_2 = c/a = 1$. Since $x_2 \approx -10^8$, the true root is $x_1 = 1 / x_2 = -10^{-8} \neq 0$.

Stabilized calculation: Compute the non-cancelling root first: $x_2 = \frac{-b - \operatorname{sgn}(b)\sqrt{b^2 - 4ac}}{2a} = -1.0 \times 10^8$ Obtain the complementary root via Vieta's formula: $x_1 = \frac{c}{a x_2} = \frac{1}{1 \cdot (-10^8)} = -1.0 \times 10^{-8}$ Plugging $x_1 = -10^{-8}$ back confirms accuracy: $(-10^{-8})^2 + 10^8(-10^{-8}) + 1 = 10^{-16} \approx 0$. The textbook formula produced a $100\%$ relative error.