Numerical Analysis

The bisection method is a simple, guaranteed-to-converge algorithm for finding a root of a continuous function f on an interval [a, b] where f(a) and f(b) have opposite signs. By the Intermediate Value Theorem, a root must lie somewhere in that bracket. The method repeatedly evaluates f at the midpoint, keeps whichever half still brackets a sign change, and discards the other half — halving the interval's width at every step. It converges linearly (each iteration gains roughly one bit of accuracy) and is slower than Newton's method, but unlike Newton's method it cannot diverge: as long as the initial bracket has a sign change, convergence is certain.

Newton's method (also called the Newton–Raphson method) is an iterative algorithm for finding a root of a differentiable function f(x) = 0. Starting from an initial guess x₀, it repeatedly replaces the function by its tangent line at the current estimate and moves to where that tangent crosses the x-axis. Near a simple root, and given a sufficiently good starting point, the method converges quadratically — the number of correct decimal digits roughly doubles with each iteration. This speed makes it the default root finder in most numerical software, though it requires the derivative f′(x) and can fail to converge if f′ is near zero or the initial guess is poor.

Fixed-point iteration solves an equation by rewriting it in the form x = g(x) and then repeatedly applying g, starting from an initial guess: x₀, x₁ = g(x₀), x₂ = g(x₁), and so on. A value x* with x* = g(x*) is called a fixed point of g, and if the iterates converge, they converge to such a fixed point — which is exactly a solution of the original equation. This is one of the simplest and most general iterative schemes in numerical analysis: Newton's method, the Babylonian square-root algorithm, and many implicit-equation solvers (like the Colebrook pipe-friction equation) are all special cases of fixed-point iteration for a particular choice of g. Whether the iteration converges, and how fast, is governed entirely by the size of g′ near the fixed point.

The condition number of a matrix A measures how much the relative error in the solution of a linear system Ax = b can be amplified by small relative errors in A or b. It is defined as κ(A) = ‖A‖·‖A⁻¹‖ for an invertible matrix and any chosen matrix norm; a condition number close to 1 indicates a well-conditioned matrix, where errors do not grow much, while a very large condition number indicates an ill-conditioned matrix, where tiny perturbations (from measurement noise or floating-point rounding) can produce wildly inaccurate solutions. Condition number is a property of the matrix itself, not of any particular algorithm — even an exact, infinite-precision solver would need extreme precision on the inputs to control the output error for an ill-conditioned system.

LU decomposition factors a square matrix A into the product of a lower triangular matrix L (with 1's on the diagonal) and an upper triangular matrix U, so that A = LU. It is the matrix formulation of Gaussian elimination: the multipliers used to eliminate entries below the diagonal during elimination become the entries of L, and the resulting row-echelon form is U. Once A is factored, solving Ax = b for many different right-hand sides b becomes cheap, since each solve reduces to a forward substitution (Ly = b) followed by a back substitution (Ux = y), each taking only O(n²) operations instead of the O(n³) needed to redo elimination from scratch. When a zero pivot would otherwise be encountered, row-permutation is added, giving PA = LU with a permutation matrix P.

Gaussian elimination is the standard systematic method for solving a system of linear equations Ax = b (or more generally for row-reducing any matrix). It transforms the augmented matrix [A | b] into upper-triangular (row echelon) form using a sequence of elementary row operations — swap two rows, scale a row by a nonzero constant, or add a multiple of one row to another — none of which change the solution set. Once the system is triangular, the last equation involves only one unknown, and back substitution unwinds the rest one variable at a time. Named after Carl Friedrich Gauss (though the technique appears far earlier, in the Chinese mathematical text The Nine Chapters on the Mathematical Art), it is the workhorse behind virtually every numerical linear-algebra routine, including LU decomposition, matrix inversion, and the structural-engineering stiffness method, and runs in O(n³) arithmetic operations for an n×n system.

Least-squares approximation finds the curve from a chosen family (most often a straight line or low-degree polynomial) that best fits a set of data points when no curve can pass through all of them exactly — which is the usual situation with noisy measurements. 'Best' means minimizing the sum of the squared vertical distances (residuals) between the data and the fitted curve, a criterion introduced by Legendre (1805) and Gauss (who used it in 1795 and later gave its probabilistic justification). For a straight-line fit this is ordinary linear regression; more generally, fitting any linear-in-its-coefficients model to overdetermined data reduces to solving a system of linear equations called the normal equations. Least squares is ubiquitous — used to fit trend lines to economic data, calibrate sensors, reconcile surveying measurements, and initialize orbit determination in astronomy.

The Jacobi and Gauss-Seidel methods are iterative alternatives to Gaussian elimination for solving Ax = b: instead of directly triangularizing the matrix, they start from a guess and repeatedly refine it until the residual is negligibly small. Both methods split A = D − L − U into its diagonal part D and its below-diagonal (−L) and above-diagonal (−U) parts, then solve each equation i for xᵢ in terms of the other current variables. Jacobi updates every component simultaneously using only values from the previous iteration, while Gauss-Seidel immediately reuses freshly updated components within the same sweep, which typically makes it converge roughly twice as fast. Both methods are only guaranteed to converge under conditions such as strict diagonal dominance of A, but when applicable — as with the large, sparse matrices that arise from discretizing partial differential equations — they are far cheaper per iteration than direct elimination, since each sweep costs only O(n²) (or O(n) for sparse systems) rather than the O(n³) of full elimination.

Finding eigenvalues exactly means finding roots of the characteristic polynomial det(A − λI) = 0, but for n larger than 4 there is no general algebraic formula for polynomial roots (a consequence of the Abel-Ruffini theorem), so every practical eigenvalue algorithm is necessarily iterative and approximate. Power iteration is the simplest such method: repeatedly multiply a vector by A and rescale, and the vector's direction converges to the eigenvector of the largest-magnitude (dominant) eigenvalue while the scaling factor converges to that eigenvalue itself. Production numerical libraries (LAPACK, and hence NumPy, MATLAB, etc.) instead use the QR algorithm, which repeatedly factors a matrix as A = QR and reassembles it as A' = RQ; this sequence of similar matrices converges to a triangular (or block-triangular) form whose diagonal entries are all the eigenvalues at once, typically after first reducing A to a cheaper-to-iterate Hessenberg form. Both families underlie everything from Google's original PageRank (power iteration on the web graph) to vibration-mode analysis in structural engineering (QR-based solvers on stiffness/mass matrices).

Simple iterative methods like Jacobi and Gauss-Seidel reduce error at every step, but they do so unevenly: they are very effective at damping high-frequency (rapidly oscillating) components of the error, yet painfully slow at removing smooth, low-frequency error components, so residual reduction stalls after a promising start on real, large-scale problems such as PDE discretizations. Multigrid methods exploit this observation directly: after a few 'smoothing' sweeps of Jacobi or Gauss-Seidel wipe out the high-frequency error, the remaining smooth error is transferred ('restricted') to a coarser grid, where it looks relatively higher-frequency and can again be smoothed away cheaply; the correction found on the coarse grid is then transferred back ('prolongated') to the fine grid. Applied recursively across a whole hierarchy of grids — the essence of the V-cycle and W-cycle schedules — multigrid can solve many elliptic PDE-derived linear systems in a number of operations that scales linearly with the number of unknowns, in stark contrast to the O(n³) of direct elimination or the slow, ever-larger iteration counts required by Jacobi/Gauss-Seidel alone on fine grids.

Euler's method is the simplest numerical technique for solving an initial value problem y′ = f(t, y), y(t₀) = y₀. It approximates the solution curve by taking small steps of size h, at each step moving in the direction given by the current slope f(t, y) as if that slope stayed constant over the whole step. Starting from the known initial point, this produces a sequence of points that traces out an approximate solution. Euler's method is simple to derive and implement but only first-order accurate (global error O(h)), so in practice it is mainly a pedagogical stepping stone to higher-order methods like Runge–Kutta, which achieve far better accuracy for the same step size.

Runge–Kutta methods are a family of iterative techniques for numerically solving initial value problems y′ = f(t, y), y(t₀) = y₀, that improve on Euler's method by sampling the slope f at several points within each step and combining them into a weighted average. The most widely used member, the classical fourth-order Runge–Kutta method (RK4), evaluates f four times per step and achieves global error O(h⁴) — dramatically more accurate than Euler's O(h) for the same step size. RK4 remains one of the most common general-purpose ODE solvers in scientific computing because of its strong accuracy-to-cost ratio and the fact that it needs no derivatives of f, only extra function evaluations.

Floating-point error is the discrepancy between an exact real-number computation and what a computer actually produces, because computers represent real numbers with only a finite number of significant binary digits (the IEEE 754 double-precision format uses about 15–17 significant decimal digits). Every arithmetic operation can introduce a small rounding error, and these errors accumulate and sometimes amplify through a long computation. The most dangerous failure mode is catastrophic cancellation: subtracting two nearly equal floating-point numbers can wipe out most of the significant digits, leaving a result dominated by rounding noise. Understanding floating-point error is essential for writing numerically stable algorithms — ones where small input or rounding errors do not blow up into large output errors.

Numerical integration, or quadrature, approximates the value of a definite integral ∫ₐᵇ f(x) dx when no closed-form antiderivative exists, or when f is known only from discrete data points. The core idea is to replace the integrand by a simple approximating function — a constant, a line, or a low-degree polynomial — over small subintervals, integrate that approximation exactly, and sum the pieces. The two most common elementary rules, the trapezoidal rule and Simpson's rule, arise from approximating f by piecewise linear and piecewise quadratic interpolants respectively. Quadrature underlies computations from areas under experimental curves to probability calculations to solving differential equations.

Simpson's rule approximates a definite integral by fitting a quadratic (parabolic) curve through each group of three consecutive sample points, rather than a straight line as in the trapezoidal rule. This extra curvature makes it exact for any polynomial up to degree 3 and gives a truncation error of O(h⁴) — far smaller than the trapezoidal rule's O(h²) for the same step size. The composite version divides [a, b] into an even number n of equal subintervals and applies alternating weights of 4 and 2 to interior nodes, with weight 1 at the endpoints. Because of its high accuracy per function evaluation, Simpson's rule is a standard default quadrature routine in numerical software.

The trapezoidal rule approximates a definite integral by replacing the curve y = f(x) with straight-line segments connecting consecutive sampled points, then summing the areas of the resulting trapezoids. It is the simplest quadrature rule beyond a single rectangle and is exact for any linear function. Applied over n equal subintervals of width h = (b−a)/n (the 'composite' trapezoidal rule), its error is O(h²): halving the step size quarters the error. Because it needs only function values (no derivatives) and is trivial to implement, it is a standard tool for integrating both formulas and tabulated experimental data.

Monte Carlo integration estimates a definite integral by averaging the integrand at random sample points rather than evaluating it on a fixed grid, as the trapezoidal or Simpson's rules do. For ∫ₐᵇ f(x) dx, the estimate is (b−a) times the average of f at N points drawn uniformly at random from [a, b]. This is a direct consequence of the law of large numbers — the sample mean of f converges to the true average of f over the interval as N grows — and the error shrinks like 1/√N regardless of the number of dimensions, unlike grid-based quadrature rules whose cost explodes exponentially with dimension (the 'curse of dimensionality'). This dimension-independence is precisely why Monte Carlo integration is the method of choice for high-dimensional integrals in computational physics, finance (option pricing), and Bayesian statistics, even though its 1/√N convergence is slower than a grid method in low dimensions.

Numerical methods are systematic recipes for computing approximate solutions to mathematical problems that have no convenient closed form — or whose exact solution is too expensive to obtain. Instead of an exact formula they produce a number, correct to a controllable tolerance, using only the four arithmetic operations a computer can perform. Root finding (solve f(x)=0), interpolation (fit a curve through data), numerical integration (quadrature), numerical differentiation, and the numerical solution of differential equations are the core tasks. They are the computational backbone of every engineering analysis package, from structural finite-element solvers to hydraulic and traffic simulations.

The finite element method (FEM) is the dominant numerical technique for solving the partial differential equations of continuum mechanics — stress and deflection in structures, heat flow, seepage, and vibration. Its central idea is discretization: a complicated body is subdivided into many small, simply-shaped pieces (the 'finite elements' — bars, triangles, quadrilaterals, tetrahedra) connected at nodes. On each element the unknown field is approximated by simple interpolation (shape) functions, the governing PDE is enforced in an averaged 'weak' sense, and each element contributes a small stiffness matrix. Assembling all of them produces one large linear system K·d = F, whose solution gives the nodal displacements (or temperatures, or heads) of the whole structure. FEM is what powers essentially every modern structural-analysis and CAE package (SAP2000, ETABS, ABAQUS, ANSYS).

Interpolation is the problem of constructing a function — almost always a polynomial — that passes exactly through a given set of data points, so that values between the known points can be estimated. Given n+1 distinct points (x₀,y₀), …, (xₙ,yₙ), there is exactly one polynomial of degree ≤ n that passes through all of them; the Lagrange form and Newton's divided-difference form are two different, algebraically equivalent recipes for writing down that same unique polynomial. Interpolation underlies reading between tabulated values (steam tables, standard normal tables), computer graphics (smooth curves through control points), and is the conceptual seed of numerical integration and differentiation, both of which approximate a function by an interpolating polynomial and then integrate or differentiate that polynomial exactly.

A single high-degree polynomial forced through many data points tends to oscillate wildly between the nodes (Runge's phenomenon), especially near the ends of the interval. Spline interpolation sidesteps this by using a different low-degree polynomial (almost always cubic) on each subinterval between consecutive data points, chosen so the pieces join up smoothly: the value, first derivative, and second derivative all match at each interior node. This piecewise construction gives a curve that looks and behaves smoothly to the eye while staying numerically well-behaved even with many data points — which is why splines are the standard tool behind font and vector-graphics curves, CAD/CAM surface design, and smooth trajectory planning in robotics and animation, as well as smooth interpolation of scientific and financial time series.

Every numerical computation carries two distinct kinds of error: input error already present in the data (measurement uncertainty, rounding from an earlier step) and the way each subsequent arithmetic operation propagates and possibly amplifies that error. A stable algorithm keeps errors from growing much faster than the problem's own inherent sensitivity (its condition number) would demand; an unstable algorithm can turn tiny, unavoidable input errors into huge output errors even when the underlying mathematical problem is perfectly well-conditioned. The most common culprit is catastrophic cancellation: subtracting two nearly equal floating-point numbers destroys most of their significant digits, since the small true difference is represented with the same absolute (not relative) precision as the large original numbers. Understanding error propagation is what separates a numerically responsible implementation of a formula from one that quietly returns garbage — the quadratic formula, variance calculations, and finite-difference derivatives are classic examples where the 'obvious' formula is numerically unstable and a reformulated version is preferred in practice.