Numerical Analysis: Algorithms for Approximate Solutions
Introduction
Numerical analysis is the study of algorithms that use numerical approximation for problems of continuous mathematics. Many mathematical problems — solving nonlinear equations, computing integrals of complicated functions, solving differential equations — cannot be solved exactly with pencil and paper. Numerical analysis provides methods that produce approximate solutions with controlled error.
The field bridges pure mathematics and computer science. Its practitioners must understand the mathematical properties of the problem, the computational properties of the algorithm, and the limitations of floating-point arithmetic on real computers. A method that works perfectly in exact arithmetic may fail catastrophically when implemented with finite precision.
Floating-Point Arithmetic and Error
Computers represent real numbers approximately using floating-point arithmetic. The IEEE 754 standard specifies formats for single and double precision, each with a finite number of bits for the significand and exponent. Most real numbers cannot be represented exactly, leading to roundoff error.
Machine Epsilon
Machine epsilon ε_m is the smallest number such that 1 + ε_m > 1 in floating-point arithmetic. For double precision (64-bit), ε_m ≈ 2.2 × 10⁻¹⁶. This value determines the relative precision of floating-point operations. Any computed result has a relative error of at least ε_m, and algorithms that amplify this error produce inaccurate results.
The condition number κ = ||A||·||A⁻¹|| of a problem measures its sensitivity to input perturbations. A problem with large condition number is ill-conditioned: small input errors produce large output errors. Even a stable algorithm cannot produce accurate results for ill-conditioned problems. Computing the condition number helps determine whether a computed result can be trusted.
Computers represent real numbers approximately using floating-point arithmetic. The IEEE 754 standard specifies formats for single and double precision, each with a finite number of bits for the significand and exponent. Most real numbers cannot be represented exactly, leading to roundoff error.
Sources of Error
Roundoff error arises from finite precision. When adding 1 and 10⁻¹⁶ in double precision, the result is exactly 1 — the small number is lost. Catastrophic cancellation occurs when subtracting nearly equal numbers, destroying significant digits. Truncation error arises from approximating infinite processes by finite ones, like replacing a derivative with a finite difference.
Stability analysis determines whether errors grow or decay as computation proceeds. A stable algorithm keeps errors under control. An unstable algorithm amplifies small errors into meaningless results. Backward error analysis relates the computed result to a slightly perturbed problem, providing insight into why some algorithms succeed and others fail.
Root Finding
Finding solutions to f(x) = 0 is a fundamental problem. The bisection method repeatedly halves an interval containing a root, guaranteeing convergence for continuous functions with opposite signs at endpoints. It converges slowly (one bit per iteration) but reliably.
Fixed-Point Iteration
Many root-finding methods can be viewed as fixed-point iterations x_{n+1} = g(x_n), where the root of f is a fixed point of g. The iteration converges if |g′(x*)| < 1 near the fixed point. Linear convergence occurs when g′(x*) ≠ 0; superlinear convergence when g′(x*) = 0.
Newton’s method is a fixed-point iteration with g(x) = x − f(x)/f′(x). It achieves quadratic convergence when f′(x*) ≠ 0 and f is sufficiently smooth. The secant method approximates f′(x) by a finite difference, achieving superlinear convergence with order approximately 1.618 without requiring derivative evaluations.
Finding solutions to f(x) = 0 is a fundamental problem. The bisection method repeatedly halves an interval containing a root, guaranteeing convergence for continuous functions with opposite signs at endpoints. It converges slowly (one bit per iteration) but reliably.
Newton’s method uses the derivative to find roots quickly: x_{n+1} = x_n − f(x_n)/f′(x_n). Near a simple root, convergence is quadratic — the number of correct digits doubles each iteration. Far from a root, Newton’s method may diverge. The secant method approximates the derivative with a finite difference, sacrificing some speed but avoiding the need for derivative calculations.
Interpolation and Approximation
Interpolation constructs a function that passes exactly through given data points. Polynomial interpolation is the classic approach: given n+1 points, there is a unique polynomial of degree at most n that passes through all of them. Lagrange and Newton forms provide numerically stable ways to compute and evaluate the interpolating polynomial.
Challenges with Polynomial Interpolation
High-degree polynomial interpolation oscillates wildly near the interval endpoints, a phenomenon observed by Runge. Fitting a single polynomial to many points typically produces poor results away from the interpolation nodes. This effect is severe for equally spaced points on functions like f(x) = 1/(1 + 25x²), where the interpolating polynomial deviates dramatically near the boundaries.
Chebyshev nodes — the roots of Chebyshev polynomials — cluster near the interval endpoints and minimize the maximum interpolation error. Using Chebyshev nodes instead of equally spaced points dramatically reduces Runge’s phenomenon. For analytic functions, Chebyshev interpolation achieves exponential convergence, making it the method of choice for high-accuracy approximation.
High-degree polynomial interpolation oscillates wildly near the interval endpoints, a phenomenon observed by Runge. Fitting a single polynomial to many points typically produces poor results away from the interpolation nodes. Spline interpolation addresses this by fitting low-degree polynomials piecewise, ensuring smoothness at the junctions.
Cubic splines are the most common choice: third-degree polynomials between each pair of points, with continuous first and second derivatives. Splines combine the flexibility of low-degree approximation with global smoothness, making them ideal for computer graphics, data fitting, and numerical differentiation.
Least Squares Approximation
When data contains noise, exact interpolation is undesirable — fitting noise does not reveal the underlying trend. Least squares approximation finds the best-fitting curve (typically a low-degree polynomial) that minimizes the sum of squared residuals. This approach underlies regression analysis and is essential for statistics and data analysis.
Numerical Integration and Differentiation
Numerical integration (quadrature) approximates definite integrals. The simplest methods approximate the integral by summing function values at equally spaced points, weighted by appropriate coefficients.
Newton-Cotes Formulas
The midpoint rule approximates the integral by the area of a rectangle at the interval midpoint. The trapezoidal rule uses a trapezoid. Simpson’s rule fits a parabola through three points, achieving higher accuracy. Romberg integration combines trapezoidal rule results at different step sizes through Richardson extrapolation.
These methods work well for smooth functions over short intervals. For oscillatory integrands or infinite intervals, specialized methods like Gaussian quadrature or adaptive quadrature are needed.
Gaussian Quadrature
Gaussian quadrature chooses evaluation points and weights to maximize accuracy for polynomials. With n points, Gaussian quadrature integrates exactly polynomials of degree up to 2n−1. The points are roots of Legendre polynomials. This method is the standard for numerical integration in scientific computing.
Numerical Differentiation
Approximating derivatives from function values is inherently ill-conditioned — subtracting nearly equal function values amplifies roundoff error. Central differences f′(x) ≈ [f(x+h) − f(x−h)]/(2h) are second-order accurate, with error proportional to h². The optimal h balances truncation error (decreasing with h) and roundoff error (increasing as 1/h).
Richardson extrapolation combines approximations at different step sizes to cancel error terms. Computing the central difference at step sizes h and h/2 and combining them gives a fourth-order accurate approximation without additional function evaluations. This technique generalizes to any method with known error expansion.
Solving Linear Systems
Gaussian elimination with partial pivoting is the standard direct method for solving Ax = b. The LU decomposition stores the elimination steps for reuse. For large sparse systems, direct methods fill in the sparse structure, making them impractical.
The condition number κ(A) = ||A||·||A⁻¹|| measures how sensitive the solution is to perturbations in A and b. A large condition number indicates an ill-conditioned system where small input errors produce large solution errors. Iterative refinement improves the solution of an ill-conditioned system by computing residuals in higher precision and correcting the solution.
Iterative methods like Jacobi, Gauss-Seidel, and conjugate gradient solve Ax = b by producing successively better approximations. Conjugate gradient methods are exceptionally efficient for large, sparse, symmetric positive-definite systems. Multigrid methods accelerate convergence by solving on multiple scales simultaneously.
Numerical Solutions of ODEs
Numerical methods for ordinary differential equations step forward in time. Euler’s method y_{n+1} = y_n + h·f(t_n, y_n) is first-order accurate. Runge-Kutta methods achieve higher order by evaluating f at intermediate points. The classic fourth-order Runge-Kutta method (RK4) is the workhorse of ODE integration.
Stiff equations — those with rapidly decaying transient components — require implicit methods that solve nonlinear equations at each step. Backward differentiation formulas (BDF) handle stiffness efficiently. Understanding stiffness is critical for chemical kinetics, electronic circuits, and other problems with multiple time scales.
What is the difference between roundoff error and truncation error? Roundoff error results from finite precision arithmetic. Truncation error results from approximating an infinite process by a finite one, such as replacing a derivative with a finite difference.
Why does Newton’s method converge so quickly? Near a simple root, Newton’s method has quadratic convergence: the error squares each iteration. Each iteration roughly doubles the number of correct digits.
When should I use spline interpolation instead of polynomial interpolation? Use splines when you have many data points to avoid polynomial oscillation. Splines provide smoothness without the wild behavior of high-degree polynomials.
What makes an ODE stiff and why does it matter? A stiff ODE has components that decay at very different rates. Explicit methods require impractically small step sizes for stiff problems. Implicit methods handle stiffness efficiently.
Computational Mathematics — Ordinary Differential Equations — Linear Algebra Applied