Numerical Analysis
Numerical analysis designs and analyses algorithms that obtain approximate solutions to mathematical problems too complicated for closed-form treatment. Every modern engineering simulation — from weather forecasts to bridge stresses to drug-protein docking — depends on it.
For a numerical method with step-size h, the error term scales as O(h^p). The integer p is the order: doubling p means each halving of h reduces the error by a factor of 2^p. Higher-order methods converge faster but may have other drawbacks (stability, complexity).
Sources of error
- Round-off error — finite-precision floating-point arithmetic. IEEE 754 double precision has ~15–16 decimal digits.
- Truncation error — discarding tail of an infinite series or higher-order terms in a Taylor expansion.
- Discretisation error — replacing continuous problems with discrete approximations (grids).
- Modelling error — equations themselves only approximate reality.
Numerical analysts seek stable, convergent and efficient algorithms — properties not guaranteed simply by being correct on paper.
Root-finding
Solving f(x) = 0 numerically. Three standard methods:
Bisection
If f is continuous and f(a)f(b) < 0, halve the interval [a, b] repeatedly. Guaranteed convergence; linear (order 1) — slow but bullet-proof.
Newton–Raphson
x_ = x_n − f(x_n)/f'(x_n)
Quadratic convergence (order 2) when f' is non-zero and the initial guess is sufficiently close. May diverge for poor starting values or near multiple roots.
Secant method
x_ = x_n − f(x_n) · (x_n − x_) / (f(x_n) − f(x_))
Superlinear order ≈ 1.618 (golden ratio); no derivative needed.
Interpolation and approximation
Given data (x_i, y_i), construct a function passing through them.
- Lagrange interpolation: p(x) = Σ y_i L_i(x) where L_i(x_j) = δ_ij. Exact at the nodes but oscillates badly (Runge phenomenon) for high degree.
- Newton's divided differences: same polynomial in a form easier to update with new points.
- Cubic splines: piecewise cubic polynomials with continuity of value, slope and curvature — smooth and stable; preferred for many applications.
- Least squares: when data are noisy, fit a polynomial of low degree minimising Σ(y_i − p(x_i))².
- For n + 1 distinct nodes there is a unique interpolating polynomial of degree ≤ n.
- High-degree polynomial interpolation on equally spaced nodes is ill-conditioned (Runge); use Chebyshev nodes or splines instead.
- Least-squares fitting minimises the L² norm of residuals; the normal equations are A^T A x = A^T b.
Numerical integration
Approximating ∫_a^b f(x) dx by sums of function values (quadrature):
| Rule | Form | Error |
|---|---|---|
| Trapezoidal | (h/2)[f(a) + 2Σf(x_i) + f(b)] | O(h²) f'' |
| Simpson's 1/3 | (h/3)[f₀ + 4f₁ + 2f₂ + … + f_n] | O(h⁴) f⁽⁴⁾ |
| Simpson's 3/8 | (3h/8)[f₀ + 3f₁ + 3f₂ + 2f₃ + …] | O(h⁴) f⁽⁴⁾ |
| Gauss–Legendre (n-pt) | Σw_i f(x_i) at optimal nodes | Exact for polys ≤ 2n − 1 |
| Romberg | Recursive trapezoidal + Richardson | High order |
Simpson's rule requires an even number of subintervals; it integrates cubics exactly thanks to symmetry.
Numerical differentiation
Forward, backward and centred finite differences:
- Forward: f'(x) ≈ [f(x+h) − f(x)]/h, O(h)
- Centred: f'(x) ≈ [f(x+h) − f(x−h)]/(2h), O(h²)
- Second derivative: f''(x) ≈ [f(x+h) − 2f(x) + f(x−h)]/h², O(h²)
Smaller h reduces truncation error but magnifies round-off — there is an optimal step size beyond which accuracy degrades.
Linear systems
For Ax = b with n equations:
- Direct methods: Gaussian elimination, LU decomposition (O(n³)), Cholesky for symmetric positive-definite A.
- Iterative methods: Jacobi, Gauss–Seidel, SOR (Successive Over-Relaxation), Conjugate Gradient.
- Conditioning: cond(A) = ‖A‖‖A⁻¹‖. Large cond → small residual ≠ small error.
The conjugate gradient method solves symmetric positive-definite systems in at most n iterations (in exact arithmetic) and is the workhorse for large sparse linear systems.
Eigenvalue problems
Av = λv. Standard algorithms:
- Power iteration finds the dominant eigenvalue.
- QR algorithm computes all eigenvalues (the de-facto standard).
- Lanczos / Arnoldi for large sparse matrices.
The QR algorithm and the Simplex method (linear programming) were ranked among the top ten algorithms of the 20th century. Other entries: FFT, Krylov subspace methods, Monte Carlo, fast multipole, Metropolis. Worth remembering for GS/CSS exams that may ask about historic algorithms.
Numerical solution of ODEs
For the IVP y' = f(t, y), y(t₀) = y₀:
| Method | Update | Order |
|---|---|---|
| Euler (explicit) | y_ = y_n + h f(t_n, y_n) | 1 |
| Improved Euler / Heun | predictor-corrector | 2 |
| RK2 (midpoint) | midpoint evaluation | 2 |
| RK4 (classical) | weighted 4 slopes | 4 |
| Adams–Bashforth | uses previous k values | k |
| Implicit / stiff methods | backward Euler, BDF, Gear | 1+ |
For stiff equations (very different timescales), implicit methods are essential despite their per-step cost.
Partial differential equations
Finite difference, finite element and finite volume methods discretise PDEs over space and time. Stability is governed by criteria such as CFL (Courant–Friedrichs–Lewy) for hyperbolic PDEs:
c Δt / Δx ≤ 1
Modern numerical PDE methods underpin computational fluid dynamics, structural analysis (FEM), seismic exploration and quantitative finance.
Modern themes
Contemporary numerical analysis emphasises:
- Sparse linear algebra for billion-variable systems.
- Iterative methods with preconditioning.
- Adaptive mesh refinement and error estimators.
- Parallel and GPU computing.
- Probabilistic numerics combining Bayesian inference with classical methods.
- Differentiable programming and scientific machine learning.
The boundary between numerical analysis, applied mathematics and computer science is becoming blurred, with the same algorithms reappearing in AI training, weather prediction and electronic-structure calculation.