Numerical Stability & Loading Calculations

1. Numerical Stability Fundamentals

Every numerical computation introduces errors. The central question of numerical stability is whether those errors grow or stay bounded as the computation proceeds. Two distinct sources of error must be understood: round-off error from finite-precision arithmetic and truncation error from replacing exact mathematical operations with finite approximations.

Round-Off Error and Floating-Point Arithmetic

IEEE 754 double-precision floating-point represents numbers as m * 2^e where m has 52 mantissa bits and e is an 11-bit exponent. Any real number that is not exactly representable is rounded to the nearest floating-point value. The machine epsilon eps_mach is approximately 2.22 × 10^-16 for double precision — the smallest number such that 1 + eps_mach is distinguishable from 1.

Truncation Error

Truncation error arises from approximating continuous mathematics with discrete or finite operations. The forward-difference approximation to a derivative is a canonical example: using f'(x) ≈ [f(x+h) - f(x)] / h truncates the Taylor series f(x+h) = f(x) + hf'(x) + (h²/2)f''(x) + ..., leaving a local truncation error of order h. Halving h halves the truncation error but doubles the number of operations, increasing accumulated round-off. The optimal h balances these two competing effects.

Condition Numbers and Ill-Conditioning

The condition number of a problem measures its inherent sensitivity to input perturbations, independent of the algorithm used. For a linear system Ax = b, the condition number kappa(A) = ||A|| ||A^-1|| bounds the relative error amplification:

An ill-conditioned problem loses approximately log10(kappa(A)) digits of accuracy in the computed solution. If kappa(A) = 10^10 and we work in double precision with ~16 significant digits, we may only get 6 correct digits in the answer — even with a perfectly stable algorithm.

2. Stability of ODE Solvers

To analyze the stability of numerical ODE methods, we use Dahlquist's test equation: dy/dt = lambda*y, where lambda is a complex number with Re(lambda) less than 0 so the true solution decays to zero. The region of absolute stability of a method is the set of step sizes h such that the numerical solution also decays to zero.

A-Stability and L-Stability

A method is A-stable if its region of absolute stability contains the entire left half-plane {Re(z) < 0} where z = h*lambda. This means the method is unconditionally stable for any decaying test equation regardless of step size. A method is L-stable if it is A-stable and additionally the amplification factor R(z) → 0 as z → -infinity. L-stability ensures that very stiff components (large negative eigenvalues) are immediately damped rather than being slowly reduced.

Stiff Equations and the Stiffness Ratio

A system y' = f(y) is stiff if the Jacobian J = df/dy has eigenvalues lambda_i with widely varying magnitudes. The stiffness ratio is max|Re(lambda_i)| / min|Re(lambda_i)|. For a stiff system with stiffness ratio 10^6, an explicit method with stability region |h*lambda| ≤ C must use h ≤ C / max|lambda| — a step size 10^6 times smaller than the accuracy requirement would dictate. This makes explicit methods prohibitively expensive for stiff problems.

Consistency and Convergence (Lax Equivalence Theorem)

For a well-posed linear initial value problem and a consistent numerical method, stability is equivalent to convergence. This is the Lax-Richtmyer equivalence theorem. Consistency means the local truncation error goes to zero as h → 0. Stability means bounded error amplification. Together they guarantee that the global error vanishes as h → 0 — the method converges to the true solution.

3. Runge-Kutta Methods

Runge-Kutta (RK) methods are one-step methods that achieve high-order accuracy by evaluating the derivative function at multiple intermediate points within each step. They are defined by a Butcher tableau that specifies the stage values and quadrature weights.

The Classical RK4 Method

The classical fourth-order Runge-Kutta method is the most widely used ODE solver in practice. It requires four function evaluations per step and achieves local truncation error of O(h^5) and global error of O(h^4).

k1 = h * f(t_n,         y_n)
k2 = h * f(t_n + h/2,   y_n + k1/2)
k3 = h * f(t_n + h/2,   y_n + k2/2)
k4 = h * f(t_n + h,     y_n + k3)

y_{n+1} = y_n + (k1 + 2k2 + 2k3 + k4) / 6

Embedded Pairs and Adaptive Step Control

Production ODE solvers use embedded Runge-Kutta pairs: two methods of different orders sharing the same function evaluations. The difference between the two solutions estimates the local error, which drives adaptive step size control. The Dormand-Prince pair (DOPRI5, used in MATLAB's ode45) embeds a 4th-order method within a 5th-order method using 6 function evaluations per step.

Implicit Runge-Kutta Methods

Explicit RK methods have limited stability regions. Implicit RK methods (like Gauss-Legendre collocation and SDIRK — Singly Diagonally Implicit Runge-Kutta) solve a system of nonlinear equations at each step but achieve A-stability or even L-stability. SDIRK methods have a diagonal Butcher tableau, so each stage can be solved sequentially with a single Newton iteration per stage.

Order Conditions and Rooted Trees

Butcher's theory of rooted trees provides a systematic way to derive the algebraic conditions that a Runge-Kutta method must satisfy to achieve a given order. For order p, all rooted trees with up to p nodes must satisfy corresponding order conditions. The number of conditions grows rapidly: 1 condition for order 1, 2 for order 2, 4 for order 3, 8 for order 4, and 17 for order 5. This is why constructing high-order explicit RK methods requires solving large systems of polynomial equations.

4. Linear Multistep Methods

Linear multistep methods (LMMs) use information from multiple previous steps to advance the solution. They can achieve high order with fewer function evaluations per step than Runge-Kutta, but require special starting procedures and are subject to the Dahlquist barrier theorem.

Adams-Bashforth Methods (Explicit)

Adams-Bashforth (AB) methods interpolate the derivative history with a polynomial and integrate it exactly to advance the solution. They are explicit (no implicit solve required) and widely used for non-stiff problems.

Adams-Moulton Methods (Implicit)

Adams-Moulton (AM) methods include the current step's derivative in the interpolation, making them implicit. They achieve one higher order than the corresponding Adams-Bashforth method with the same number of prior steps and have larger stability regions. In practice, AB and AM methods are paired as Predictor-Corrector (PECE) schemes.

BDF Methods for Stiff ODEs

Backward Differentiation Formulas (BDF) are the method of choice for stiff ODEs. Rather than interpolating the derivative, BDF methods differentiate the interpolating polynomial of the solution history and set it equal to the derivative. BDF1 through BDF6 are zero-stable and A(alpha)-stable for orders 1–6.

5. Matrix Stability: Eigenvalue Analysis

The stability of iterative algorithms, linear recurrences, and dynamical systems depends critically on the eigenvalue structure of the associated matrices. The spectral radius and Gershgorin's theorem provide algebraic tools for bounding eigenvalues without full diagonalization.

Spectral Radius and Matrix Norms

The spectral radius of a matrix A is rho(A) = max|lambda_i|, the largest absolute value of any eigenvalue. For a linear recurrence x_{k+1} = Ax_k, the sequence converges to zero if and only if rho(A) less than 1. The spectral radius bounds induced matrix norms: rho(A) ≤ ||A|| for any induced norm. For symmetric matrices, rho(A) = ||A||_2 (the spectral norm equals the largest singular value).

The Gershgorin Circle Theorem

The Gershgorin circle theorem localizes eigenvalues using only the matrix entries. For a matrix A, define the ith Gershgorin disk as: D_i = { z in C : |z - A_{ii}| ≤ R_i } where R_i is the sum of absolute values of off-diagonal entries in row i. Every eigenvalue of A lies in the union of these disks.

A = [ 4   1  -1 ]
    [ 0   3   0.5]
    [-0.5  0   2 ]

D1: center 4,  radius 1 + 1 = 2  → eigenvalue in [2, 6]
D2: center 3,  radius 0 + 0.5 = 0.5 → eigenvalue in [2.5, 3.5]
D3: center 2,  radius 0.5 + 0 = 0.5 → eigenvalue in [1.5, 2.5]

All disks in right half-plane → matrix is stable (negative of A stable)

Lyapunov Stability for Continuous Systems

For the continuous system dx/dt = Ax, the origin is asymptotically stable if and only if all eigenvalues of A have strictly negative real parts (A is Hurwitz stable). This can be tested without computing eigenvalues using the Routh-Hurwitz criterion for the characteristic polynomial, or by solving the Lyapunov equation A^T P + PA = -Q for a positive definite P, given any positive definite Q.

6. LU Decomposition with Pivoting

Gaussian elimination factorizes a matrix A as A = LU where L is unit lower triangular and U is upper triangular. This factorization allows efficient solution of multiple right-hand sides by forward and backward substitution. Pivoting strategies are essential for numerical stability.

Partial Pivoting

Partial pivoting interchanges rows before each elimination step so that the pivot element has the largest absolute value in its column. This bounds all multipliers: |L_{ij}| ≤ 1 for all i > j. The factorization produces PA = LU where P is a permutation matrix.

for k = 1 to n-1:
  p = argmax_{i>=k} |A[i,k]|   // find pivot row
  swap rows k and p in A        // row permutation
  P[k] = p                      // record permutation

  for i = k+1 to n:
    L[i,k] = A[i,k] / A[k,k]   // compute multiplier
    A[i,k:n] -= L[i,k]*A[k,k:n] // eliminate

Growth Factor and Stability Guarantees

The growth factor g_n = max|U_{ij}| / max|A_{ij}| measures how much the matrix entries grow during elimination. Large growth factors indicate potential instability. With partial pivoting, the theoretical bound is g_n ≤ 2^(n-1), but in practice the growth factor is rarely larger than n. With complete pivoting (permuting both rows and columns), the bound is g_n ≤ (n * 2 * 3^(1/2) * ... * n^(1/(n-1)))^(1/2), which grows much more slowly.

Cholesky Decomposition for Symmetric Positive Definite Matrices

When A is symmetric positive definite (SPD), the Cholesky decomposition A = LL^T is always stable without pivoting, and costs half as much as LU (roughly n^3/6 flops vs n^3/3 flops). The growth factor is exactly 1: max|L_{ij}|^2 ≤ max|A_{ii}|, so entries never grow during factorization. Cholesky is the preferred method for SPD systems arising in finite element analysis, statistics (normal equations), and optimization.

for j = 1 to n:
  L[j,j] = sqrt(A[j,j] - sum(L[j,k]^2, k=1..j-1))
  for i = j+1 to n:
    L[i,j] = (A[i,j] - sum(L[i,k]*L[j,k], k=1..j-1)) / L[j,j]

// Solve Ax = b via forward sub (Ly = b) then back sub (L^T x = y)

7. Classical Iterative Solvers

For large sparse systems where direct factorization is too expensive due to fill-in, iterative methods build a sequence of approximations that converge to the solution. The simplest are the splitting methods: Jacobi, Gauss-Seidel, and SOR.

Matrix Splitting Framework

A matrix splitting writes A = M - N where M is easy to invert. The iterative scheme is Mx_{k+1} = Nx_k + b, or equivalently x_{k+1} = M^-1Nx_k + M^-1b. The iteration converges if and only if rho(M^-1N) less than 1. The asymptotic convergence rate is -log(rho(M^-1N)) per iteration.

Jacobi Method

In the Jacobi method, each component of x is updated using only values from the previous iteration. The update for component i is: x_i^{(k+1)} = (b_i - sum_{j ≠ i} A_{ij} x_j^{(k)}) / A_{ii}. Because all components are updated simultaneously, Jacobi is easily parallelized. It converges if A is strictly diagonally dominant: |A_{ii}| > sum_{j≠i} |A_{ij}|.

Gauss-Seidel Method

Gauss-Seidel uses the most recently computed component values immediately: x_i^{(k+1)} = (b_i - sum_{j<i} A_{ij} x_j^{(k+1)} - sum_{j>i} A_{ij} x_j^{(k)}) / A_{ii}. This typically converges roughly twice as fast as Jacobi for the same problem class. For symmetric positive definite systems, Gauss-Seidel always converges.

SOR and Optimal Relaxation

Successive Over-Relaxation (SOR) accelerates Gauss-Seidel by combining the current Gauss-Seidel update with the previous iterate: x_i^{new} = (1-omega)*x_i^{old} + omega*x_i^{GS}. For 0 < omega < 1 (under-relaxation) convergence is guaranteed for more problems; for 1 < omega < 2 (over-relaxation) convergence is accelerated. For the model problem (Poisson equation on a grid), the optimal omega is known analytically.

8. Krylov Subspace Methods

Krylov subspace methods are the state of the art for large sparse linear systems. They build an approximation to the solution from the Krylov subspace K_k(A, r_0) = span{r_0, Ar_0, A^2r_0, ..., A^{k-1}r_0}where r_0 = b - Ax_0 is the initial residual. Each iteration requires only one matrix-vector product, making them feasible for matrices with millions of unknowns.

Conjugate Gradient (CG) for Symmetric Positive Definite Systems

The Conjugate Gradient method is optimal for SPD systems: it minimizes the A-norm of the error ||e_k||_A = sqrt(e_k^T A e_k) over all vectors in x_0 + K_k. CG converges in at most n iterations (in exact arithmetic), but in practice converges much faster due to eigenvalue clustering.

r_0 = b - A*x_0,  p_0 = r_0
for k = 0, 1, 2, ...:
  alpha_k = (r_k^T r_k) / (p_k^T A p_k)
  x_{k+1} = x_k + alpha_k * p_k
  r_{k+1} = r_k - alpha_k * A*p_k
  beta_k  = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)
  p_{k+1} = r_{k+1} + beta_k * p_k

GMRES for Non-Symmetric Systems

GMRES (Generalized Minimal Residual) handles non-symmetric systems by minimizing ||b - Ax_k||_2 over x_0 + K_k. It builds an orthonormal basis for K_k using the Arnoldi process and solves a small least-squares problem at each step. GMRES is guaranteed to converge (the residual is monotonically non-increasing) but may require restarts to bound memory use.

Preconditioning Strategies

A preconditioner M approximates A^-1 cheaply. Left preconditioning solves M^-1Ax = M^-1b; right preconditioning solves AM^-1(Mx) = b. The goal is to cluster the eigenvalues of M^-1A (or AM^-1) near 1, dramatically reducing the number of Krylov iterations required.

9. Stiff Systems in Detail

Stiff ODEs appear throughout science and engineering wherever multiple processes operate on widely different time scales. Solving them efficiently requires specialized implicit methods and careful linearization strategies.

Stiffness Detection and Diagnosis

A system y' = f(y) is stiff near a solution if the Jacobian J = df/dy has eigenvalues lambda_i satisfying: some |Re(lambda_i)| is large while the solution of interest varies slowly. The stiffness ratio S = max|Re(lambda_i)| / min|Re(lambda_i)| quantifies severity. Practically, stiffness is detected when an explicit solver takes many more steps than accuracy would require.

Linearization and Newton Iteration for Implicit Methods

Implicit methods require solving a nonlinear system G(Y) = 0 at each step. Newton's method is used: Y^(k+1) = Y^(k) - [G'(Y^(k))]^(-1) G(Y^(k)). For BDF methods, G'(Y) = I - h*beta_s*J where J is the Jacobian of f. The Jacobian is typically computed once at the start of the step and reused for several steps (modified Newton or Jacobian lagging).

// BDF2: (3y_{n+1} - 4y_n + y_{n-1})/2h = f(y_{n+1})
// Rewrite as G(y_{n+1}) = 0:
G(y) = y - (2h/3)*f(y) - (4y_n - y_{n-1})/3

// Newton iteration:
J_G = I - (2h/3) * J_f
for k = 0, 1, ...:
  J_G * delta = -G(y^(k))   // solve linear system
  y^(k+1) = y^(k) + delta
  if ||delta|| < tol: break

Differential-Algebraic Equations (DAEs)

DAEs combine differential equations with algebraic constraints: F(t, y, y') = 0. The index of a DAE measures how many differentiations are needed to convert it to a pure ODE. Index-1 DAEs (like those from circuit simulation or incompressible Navier-Stokes) can be handled by BDF methods. Higher-index DAEs require reformulation or specialized solvers to avoid index reduction instabilities.

10. Error Analysis: Forward, Backward, and Mixed

Numerical analysts distinguish three fundamental types of error analysis. Each gives different insight into algorithm behavior and is appropriate in different contexts.

Forward Error Analysis

Forward error analysis tracks errors as they propagate through each arithmetic operation. Starting from input errors (round-off in representing data), it bounds the final output error by accumulating error bounds at each step. The forward error bound is ||computed_x - true_x|| / ||true_x||. Forward analysis tends to produce pessimistic (overly large) bounds because it assumes worst-case error accumulation at every step.

Backward Error Analysis (Wilkinson's Approach)

Backward error analysis asks: what perturbation to the original problem would produce the computed answer exactly? For Gaussian elimination with partial pivoting, Wilkinson showed that the computed solution x_computed exactly solves (A + delta_A) x_computed = b where ||delta_A||_inf / ||A||_inf ≤ g_n * eps_mach * n. This is a small backward error when the growth factor g_n is small — meaning the algorithm is backward stable even if the forward error could be large.

Mixed Stability and the Error Relationship

The fundamental relationship connecting backward error, forward error, and conditioning is:

Global vs. Local Truncation Error for ODEs

For ODE methods, the local truncation error (LTE) measures the defect introduced in a single step assuming exact initial data. An s-stage Runge-Kutta method of order p has LTE = O(h^{p+1}). The global truncation error (GTE) accumulates LTE over n = T/h steps. Since each step introduces O(h^{p+1}) error and there are O(1/h) steps, the GTE is O(h^p). This is why order-p methods achieve O(h^p) global accuracy.

11. Applications in Science and Engineering

Numerical stability analysis is not purely theoretical — it determines whether large-scale simulations give trustworthy results. The following applications illustrate the real-world stakes of choosing stable algorithms.

12. Practice Problems with Solutions

Work through these problems to consolidate your understanding of numerical stability. Expand each solution only after attempting the problem.

Quick Reference: Methods and Stability Properties