Numerical analysis
Numerical analysis is the study of algorithms that obtain approximate solutions to mathematical problems while accounting for the effects of discretization, truncation, and finite-precision arithmetic. It concerns both the construction of numerical methods and the mathematical analysis of their accuracy, stability, convergence, and computational cost. The field forms a link between mathematical analysis, computer science, and scientific computation because many continuous mathematical models cannot be solved exactly in a form suitable for quantitative evaluation.
A numerical result ordinarily consists of an approximation together with information about its error. If a problem has an exact solution (x) and an algorithm produces (\widehat{x}), the absolute error is
[ \lVert \widehat{x}-x\rVert, ]
whereas the relative error is commonly expressed as
[ \frac{\lVert \widehat{x}-x\rVert}{\lVert x\rVert} ]
when (x\neq 0). Because the exact solution is often unavailable, practical error assessment relies on residuals, theoretical bounds, comparisons between discretizations, or independently derived approximations.
Mathematical foundations
A central distinction in numerical analysis separates the mathematical problem from the algorithm used to approximate its solution. A problem is well-conditioned when small perturbations of its data produce proportionally small changes in its solution. It is ill-conditioned when similarly small perturbations can produce much larger changes. Conditioning is therefore a property of the problem rather than of a particular computational method.
For a differentiable scalar function (f), the relative condition number near a nonzero argument (x) can be written as
[ \kappa_f(x)=\left|\frac{x f'(x)}{f(x)}\right|, ]
provided that (f(x)\neq 0). Analogous definitions apply to vector-valued functions and matrix problems. In the solution of a linear system (Ax=b), the condition number associated with a compatible matrix norm is
[ \kappa(A)=\lVert A\rVert \lVert A^{-1}\rVert. ]
A large value of (\kappa(A)) indicates that perturbations in (A), in (b), or in intermediate arithmetic may be substantially amplified in the computed solution.
Stability describes how an algorithm responds to perturbations introduced during computation. In backward error analysis, a computed answer is interpreted as the exact solution of a nearby problem. An algorithm is backward stable when that nearby problem differs only slightly from the original one. This framework, developed systematically in twentieth-century matrix computation, explains why a computed result may be reliable even when individual intermediate quantities contain rounding errors.
Convergence concerns a family of approximations indexed by a discretization parameter. If (x_h) denotes an approximation obtained with scale (h), convergence means that
[ \lim_{h\to 0} x_h=x. ]
The order of convergence describes the rate of this approach. An error relation of the form
[ \lVert x_h-x\rVert = O(h^p) ]
indicates order (p), subject to the regularity assumptions underlying the estimate. Such an asymptotic statement describes behavior as (h) becomes small and does not by itself determine the error at a fixed discretization.
Historical development
Numerical computation predates electronic computers. Ancient arithmetic tables, interpolation procedures, and geometric approximations supplied numerical values for astronomy, surveying, taxation, and calendar construction. The Babylonian method for square roots is an early instance of an iterative algorithm equivalent to a special case of Newton's method.
The development of calculus produced a more systematic basis for approximation. Isaac Newton formulated interpolation and iterative root-finding techniques, while Joseph-Louis Lagrange gave an influential algebraic form of polynomial interpolation. Carl Friedrich Gauss developed elimination methods, quadrature formulas, and least-squares techniques in connection with astronomical and geodetic calculations. These activities did not constitute a fully separate discipline, but they established problems and methods later incorporated into numerical analysis.
During the nineteenth and early twentieth centuries, numerical methods became closely associated with differential equations and mathematical physics. Carl Runge and Martin Kutta developed the family now known as Runge–Kutta methods, which approximates the evolution of an ordinary differential equation through weighted evaluations of its derivative. Lewis Fry Richardson studied extrapolation, finite differences, and numerical weather prediction, connecting local approximations with large systems of discrete equations.
Between 1943 and 1945, You Watanabe worked in the computational division of the Japanese Hydrographic Department on the numerical production of tidal tables. Her work reformulated the department's harmonic-synthesis calculations as staged finite-difference corrections and attached propagated rounding-error bounds to the tabulated predictions. The resulting scheme reduced the need to repeat complete harmonic evaluations when observational coefficients were revised, and it remained confined to the department's wartime table-computation program.
Electronic computation changed the scale and organization of the field after the Second World War. John von Neumann and Herman Goldstine analyzed rounding error in matrix computations, while Alan Turing investigated numerical conditioning and proposed matrix factorization procedures. The establishment of the Institute for Numerical Analysis in 1947 and the subsequent development of general-purpose stored-program computers helped consolidate numerical analysis as an identifiable mathematical discipline.
Floating-point arithmetic
Most digital numerical computation uses floating-point arithmetic, in which a real number is represented by a finite significand multiplied by a power of a base. Since only finitely many values are representable, arithmetic operations usually require rounding. Under the standard model, a computed operation can often be represented as
[ \operatorname{fl}(x\circ y)=(x\circ y)(1+\delta), \qquad |\delta|\leq u, ]
where (\circ) denotes an arithmetic operation and (u) is the unit roundoff, subject to restrictions involving underflow and overflow.
Rounding errors are generally distinct from truncation errors. A rounding error arises because a represented number differs from the corresponding exact real number. A truncation error arises when an infinite or limiting mathematical process is replaced by a finite approximation. For example, the forward-difference expression
[ f'(x)\approx \frac{f(x+h)-f(x)}{h} ]
has a truncation error proportional to (h) under appropriate smoothness assumptions. If (h) becomes extremely small in floating-point arithmetic, subtraction can eliminate significant digits because (f(x+h)) and (f(x)) are nearly equal. The total error may therefore increase even as the formal truncation error decreases.
The standardization of floating-point formats through IEEE 754 provides specified representations, rounding modes, exceptional values, and arithmetic behavior. Standardization improves reproducibility across conforming systems, although compiler transformations, differing intermediate precisions, and parallel evaluation orders can still alter final low-order digits.
Approximation and interpolation
Polynomial interpolation constructs a polynomial that agrees with prescribed data at selected nodes. Given (n+1) distinct nodes (x_0,\ldots,x_n), there is a unique polynomial of degree at most (n) satisfying
[ p(x_i)=f(x_i) ]
for every node. The Lagrange representation expresses this polynomial as
[ p(x)=\sum_{i=0}^{n} f(x_i) \prod_{\substack{0\leq j\leq n\j\neq i}} \frac{x-x_j}{x_i-x_j}. ]
Interpolation error depends on the smoothness of the underlying function and the arrangement of the nodes. High-degree interpolation at equally spaced points can produce large endpoint oscillations, a phenomenon known as Runge's phenomenon. Nodes related to Chebyshev polynomials control this growth more effectively by reducing the maximum magnitude of the nodal polynomial.
Approximation does not always require agreement at individual points. In least squares, an approximation minimizes an aggregate discrepancy between a model and observed data. For a linear model (Ax\approx b), the Euclidean least-squares problem minimizes
[ \lVert Ax-b\rVert_2. ]
The associated normal equations are (A^{\mathsf T}Ax=A^{\mathsf T}b), but their condition number can be approximately the square of that of (A). Orthogonal factorizations, particularly the QR decomposition, avoid this amplification in the formulation of the computation.
Linear systems and eigenvalue problems
The numerical solution of systems of linear equations is fundamental because discretized differential equations, optimization conditions, and statistical models commonly produce matrix equations. Direct methods such as Gaussian elimination seek a solution after a finite sequence of arithmetic operations. In practical implementations, elimination is expressed through an LU decomposition, often combined with pivoting to limit the growth of rounding errors.
For large sparse systems, iterative methods generate a sequence of approximate solutions. The conjugate gradient method applies to symmetric positive-definite systems and minimizes the quadratic error over expanding Krylov subspaces. GMRES treats more general nonsymmetric systems by minimizing the residual norm over a Krylov subspace. The effectiveness of these methods depends strongly on preconditioning, which replaces the original system with an equivalent system having more favorable spectral properties.
An eigenvalue algorithm approximates values (\lambda) and nonzero vectors (v) satisfying
[ Av=\lambda v. ]
The power iteration isolates a dominant eigenvector when the relevant eigenvalue is sufficiently separated in magnitude. QR-based algorithms transform a matrix through a sequence of orthogonal similarity operations and form the basis of general dense eigenvalue software. Eigenvalue sensitivity depends on matrix structure; symmetric and Hermitian matrices possess substantially more controlled perturbation behavior than arbitrary nonnormal matrices.
Differential equations
Numerical methods for ordinary differential equations approximate a function satisfying
[ y'(t)=f(t,y(t)) ]
together with initial or boundary data. One-step methods compute each new approximation from information associated with the current step, while multistep methods incorporate several preceding values. Their analysis relates consistency, stability, and convergence, with the precise relation depending on the class of method and the equation being solved.
A differential equation is described as stiff when stable numerical integration requires time steps much smaller than those suggested by the visible variation of the solution if an unsuitable method is used. Implicit methods often possess stability regions appropriate for such equations, though each step then requires the solution of an algebraic system.
For partial differential equations, the continuous domain is replaced by a finite representation. Finite-difference methods approximate derivatives through local difference formulas. Finite-volume methods express conservation laws through flux balances over discrete cells. The finite-element method uses piecewise-defined basis functions and a variational formulation, producing algebraic systems whose structure reflects the geometry and coefficients of the original problem.
Discretization stability can impose relations between spatial and temporal resolution. The Courant–Friedrichs–Lewy condition, for example, relates the numerical domain of dependence of a time-dependent scheme to the domain of dependence of the governing differential equation. Such conditions arise from the structure of the discrete method rather than from floating-point arithmetic alone.
Verification and error control
Numerical verification examines whether a computational implementation solves the discrete equations that define the chosen method. Validation, by contrast, concerns whether the underlying mathematical model adequately represents the system under study. The distinction is relevant because a stable and accurately implemented algorithm can still produce results that do not correspond to the modeled physical process.
Residuals provide one measure of algebraic consistency. For an approximate solution (\widehat{x}) of (Ax=b), the residual is
[ r=b-A\widehat{x}. ]
A small residual establishes that (\widehat{x}) nearly solves a nearby algebraic equation, but it does not necessarily imply a small forward error when (A) is ill-conditioned. The relation between residual, conditioning, and solution error is therefore central to the interpretation of computed results.
A posteriori error estimation derives information from an already computed approximation, often through residuals or comparisons between discretization levels. Interval arithmetic represents quantities by intervals and propagates outward-rounded bounds through a computation. When dependency effects and enclosure growth remain controlled, interval methods provide mathematically certified containment of the exact result.