Numerical methods form the computational discipline dedicated to approximating solutions to continuous mathematical problems using finite, discrete computing machines. Because digital computers represent real numbers with finite bit lengths, exact mathematical solutions to differential equations, large linear systems, and continuous optimization problems are rarely achievable.
Numerical analysis provides the mathematical framework for constructing algorithms that balance computational complexity, numerical stability, and rigorous truncation error bounds.
In IEEE 754 double-precision arithmetic (binary64), a real number x is stored across 64 bits: 1 sign bit (s), 11 exponent bits (e), and 52 fraction/mantissa bits (m):
The machine epsilon \epsilon_{\text{mach}} = 2^{-52} \approx 2.2204 \times 10^{-16} represents the upper bound on the relative error due to rounding:
IEEE 754 Double-Precision Layout (64-bit):
+---+----------------------+----------------------------------------------------+
| s | Exponent (11 bits) | Significand / Mantissa (52 bits) |
+---+----------------------+----------------------------------------------------+
63 62 52 51 0
For a continuous mathematical problem f(x) evaluated via a numerical algorithm \hat{f}(x):
The Condition Number \kappa measures the intrinsic sensitivity of the mathematical problem to perturbations, independent of the algorithm:
For a linear system A x = b, the condition number with respect to matrix inversion is:
where \sigma_{\max}, \sigma_{\min} are the maximum and minimum singular values of A. A rule of thumb is that if \kappa(A) \approx 10^k, solving A x = b in double precision loses approximately k decimal digits of accuracy.
Catastrophic cancellation occurs when subtracting two nearly equal floating-point numbers x \approx y, amplifying the relative error of the lower-order bits. For example, evaluating roots of a x^2 + b x + c = 0 when b^2 \gg 4ac:
Computing x_2 = \frac{-b + \sqrt{b^2-4ac}}{2a} directly induces severe cancellation, whereas the algebraic reformulations above preserve full machine precision.
Solving large systems of linear equations A x = b is the computational bottleneck in finite element analysis (FEA), machine learning backpropagation, and fluid dynamics simulations.
Direct vs. Iterative Solvers:
--------------------------------------------------------------------------------
System Type Recommended Method Complexity
--------------------------------------------------------------------------------
Dense (Small/Medium) LU Factorization (Pivoting) O(2/3 n³)
Symmetric Pos. Definite Cholesky Factorization (L L^T) O(1/3 n³) [Twice as fast as LU]
Overdetermined (Least Sq)QR via Householder Reflections O(2 m n² - 2/3 n³)
Large Sparse Symmetric Preconditioned Conjugate Gradient O(k · nnz(A))
Large Sparse Asymmetric GMRES / BiCGSTAB O(k · nnz(A) + k² n)
--------------------------------------------------------------------------------
For large sparse matrices (n > 10^6), direct \mathcal{O}(n^3) factorizations require excessive memory. Krylov subspace methods project A onto the subspace \mathcal{K}_k(A, r_0) = \operatorname{span}\{r_0, A r_0, A^2 r_0, \dots, A^{k-1} r_0\}.
The Conjugate Gradient (CG) algorithm minimizes the energy norm \|x_k - x^*\|_A over \mathcal{K}_k:
import numpy as np
def conjugate_gradient(A: np.ndarray, b: np.ndarray, x0: np.ndarray = None, tol: float = 1e-10, max_iter: int = 1000):
"""
Solves A x = b for Symmetric Positive Definite (SPD) matrix A using Conjugate Gradient.
"""
n = len(b)
x = np.zeros(n) if x0 is None else x0.copy()
r = b - A @ x
p = r.copy()
rs_old = np.dot(r, r)
for i in range(max_iter):
Ap = A @ p
alpha = rs_old / np.dot(p, Ap)
x += alpha * p
r -= alpha * Ap
rs_new = np.dot(r, r)
if np.sqrt(rs_new) < tol:
return x, i + 1
p = r + (rs_new / rs_old) * p
rs_old = rs_new
return x, max_iter
Convergence rate is governed by the condition number \kappa(A):
Preconditioning transforms the system into M^{-1} A x = M^{-1} b such that \kappa(M^{-1} A) \approx 1.
To solve f(x) = 0 for f: \mathbb{R} \to \mathbb{R}:
Theorem (Quadratic Convergence): If f \in C^2, f'(x^*) \neq 0, and x_0 is sufficiently close to x^*, the error e_k = x_k - x^* satisfies:
This doubles the number of correct significant figures at every iteration near the root.
Convergence Comparison for Root Finding:
--------------------------------------------------------------------------------
Method Order of Convergence (p) Function Evals / Step
--------------------------------------------------------------------------------
Bisection p = 1 (Linear) 1 eval (Guaranteed robust)
Secant Method p = (1 + √5)/2 ≈ 1.618 1 eval (No derivative needed)
Newton-Raphson p = 2 (Quadratic) 2 evals (f and f')
Halley's Method p = 3 (Cubic) 3 evals (f, f', and f'')
--------------------------------------------------------------------------------
For unconstrained minimization \min_{x \in \mathbb{R}^n} f(x), pure Newton updates x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k) require \mathcal{O}(n^3) operations per step to compute and invert the Hessian \nabla^2 f(x).
The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method iteratively updates an inverse Hessian approximation H_k \approx (\nabla^2 f(x_k))^{-1} with rank-2 matrix updates:
This maintains positive definiteness (H_{k+1} \succ 0) and guarantees superlinear convergence with only \mathcal{O}(n^2) complexity per step.
Numerical integration approximates definite integrals I(f) = \int_a^b f(x) dx \approx \sum_{i=0}^n w_i f(x_i).
Quadrature Rules:
+-------------------+-----------------------------------+-----------------------+
| Method | Formula / Stencil | Error Order |
+-------------------+-----------------------------------+-----------------------+
| Trapezoidal Rule | h/2 [f(a) + 2∑f(x_i) + f(b)] | O(h²) · f''(ξ) |
| Simpson's 1/3 Rule| h/3 [f_0 + 4∑f_odd + 2∑f_even + fn]| O(h⁴) · f⁽⁴⁾(ξ) |
| Gauss-Legendre | ∑ w_i f(x_i) on Legendre roots | Exact for poly ≤ 2n-1 |
+-------------------+-----------------------------------+-----------------------+
By choosing non-equispaced evaluation nodes x_i as the roots of the n-th degree Legendre polynomial P_n(x) on [-1, 1], Gaussian quadrature achieves the maximum possible algebraic degree of exactness: an n-point rule integrates all polynomials up to degree 2n - 1 exactly.
Solving the initial value problem (IVP):
The general s-stage explicit Runge-Kutta method computes:
For the classical 4th-order Runge-Kutta (RK4):
RK4 has a local truncation error of \mathcal{O}(h^5) and global error of \mathcal{O}(h^4).
def rk4_step(f, t: float, y: np.ndarray, h: float) -> np.ndarray:
"""Executes a single 4th-order Runge-Kutta step."""
k1 = f(t, y)
k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
k4 = f(t + h, y + h * k3)
return y + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
A system of ODEs y' = A y is stiff if the eigenvalues \lambda_i of the Jacobian J = \frac{\partial f}{\partial y} have negative real parts with vastly different magnitudes:
Applying explicit forward Euler y_{n+1} = (1 + h \lambda) y_n to the test equation y' = \lambda y (\operatorname{Re}(\lambda) < 0) requires |1 + h \lambda| \le 1 for numerical stability, forcing an impractically small step size:
Stability Regions in Complex z = hλ Plane:
Explicit Forward Euler Implicit Backward Euler
(Bounded Circle) (Unbounded Exterior Region)
| Im(z) | Im(z)
| |
+---+---+ XXXX|XXXX
/ | \ XXXX|XXXX
| Stable | Re(z) XXXX+---+XXXX Re(z)
\ Region / XXXX| |XXXX (Stable everywhere
+---+---+ XXXX|XXXX except inside circle)
| |
An ODE solver is A-stable if its stability region includes the entire left half of the complex plane \mathbb{C}^- = \{z \in \mathbb{C} : \operatorname{Re}(z) \le 0\}.
[ Numerical Problem ]
|
+------------------+--------------+------------------+------------------+
| | | |
[Linear System] [Nonlinear Opt] [Integration] [ODEs]
| | | |
Sparse / SPD? Hessian cheap? Smooth bounds? Stiff system?
/ \ / \ / \ / \
Yes No Yes No Yes No Yes No
/ \ / \ / \ / \
[PCG] [LU / [Newton] [BFGS / [Gauss- [Adaptive [Implicit [Adaptive
GMRES] L-BFGS] Legendre] Simpson] BDF/Radau] RK45]