Numerical analysis - Major Computational Techniques

Numerical Analysis: Computing, Solving, and Optimizing Introduction Numerical analysis is the study of algorithms for solving mathematical problems using numerical approximations rather than exact analytical solutions. Many real-world problems cannot be solved by hand or do not have closed-form solutions, making numerical methods essential tools. This guide covers the major problem classes that numerical analysis addresses: evaluating functions, solving equations, analyzing matrix properties, optimizing outcomes, integrating functions, and simulating differential equations. Computing Values of Functions One of the most practical problems in numerical analysis is determining function values when you don't have a complete picture of the function itself. Often you have only discrete measurements or data points and need to find values in between them or beyond what you know. Interpolation Interpolation estimates function values at points that lie between your known data points. Imagine you have temperature measurements taken at noon and at 2 PM, but you need to estimate the temperature at 1 PM. Linear interpolation is the simplest approach: you assume the temperature changed at a constant rate between the two measurements and draw a straight line connecting them. More generally, interpolation becomes crucial when dealing with expensive or difficult measurements. If you have measurements at specific points, you can construct a smooth curve passing through all of them to estimate values elsewhere. The most common methods create polynomial curves through your data points. Extrapolation Extrapolation extends the same principle beyond your known data. Instead of predicting between measured points, you predict outside the measured range. A government agency might use historical GDP data to project future economic growth, or a physicist might use laboratory measurements to infer what happens in extreme conditions. Extrapolation is generally less reliable than interpolation because you're venturing into unknown territory without the stability of surrounding data. Regression and Least-Squares Fitting Real-world measurements are noisy and imperfect. When your data contains measurement errors or natural variation, you don't want a curve that passes exactly through every noisy point—that would amplify the noise. Instead, regression fits a smooth model (often a low-degree polynomial or linear function) to your data. The least-squares method is the standard approach. It finds the model that minimizes the sum of squared differences between the predicted values and the actual measurements. Squaring the differences ensures that both large positive and negative errors are penalized equally, and large errors are penalized more heavily than small ones. This method produces the "best fit" curve in a well-defined mathematical sense. Polynomial Evaluation: The Horner Scheme Once you have a polynomial to evaluate—say $p(x) = a0 + a1 x + a2 x^2 + a3 x^3$—you might compute it directly by calculating each power of $x$ separately. But this is computationally wasteful. The Horner scheme reorganizes the calculation as: $$p(x) = a0 + x(a1 + x(a2 + x(a3)))$$ This nested form requires far fewer multiplications and additions. For a polynomial of degree $n$, Horner's method needs only $n$ multiplications instead of $1 + 2 + \ldots + n$. When evaluating polynomials thousands or millions of times, this efficiency matters significantly. Optimization Optimization addresses the problem: "What input value produces the best output?" This could mean maximizing profit given production constraints, minimizing material usage, or finding the operating conditions that produce the highest efficiency. You're searching through possible input values to find the one that optimizes your objective. Different problems require different approaches. Simple unimodal problems (where there's a single peak) can be solved with straightforward search methods. More complex problems with multiple constraints or local optima require sophisticated algorithms. Solving Equations and Systems of Equations A fundamental problem across science and engineering is solving equations: finding the value of $x$ that satisfies $f(x) = 0$, or finding the values $(x1, x2, \ldots, xn)$ that simultaneously satisfy multiple equations. Linear systems are common and have well-developed solutions, while nonlinear equations require different approaches. Solving Linear Systems: Direct Methods A linear system has the form $Ax = b$, where $A$ is a known matrix, $b$ is a known vector, and $x$ is the unknown solution vector you seek. This appears constantly in applications. Gaussian elimination is the classical method taught in algebra. You perform row operations (swapping rows, adding multiples of one row to another) to transform the system into an upper triangular form, then solve using back-substitution. It's reliable, widely understood, and works well for moderately-sized systems. More specialized decomposition methods are efficient for certain matrix structures: LU decomposition factors $A$ into lower and upper triangular matrices ($L$ and $U$). Once computed, this decomposition lets you solve multiple systems $Ax = b1$, $Ax = b2$, etc., with much less work than repeating Gaussian elimination. Cholesky decomposition applies specifically to symmetric positive-definite matrices (matrices that equal their own transpose and satisfy certain mathematical properties). It exploits this special structure to be faster and more numerically stable than general methods. QR decomposition factors $A$ into an orthogonal matrix $Q$ and an upper triangular matrix $R$. This method handles non-square matrices and is more numerically stable than Gaussian elimination for ill-conditioned problems (problems where small changes in the input cause large changes in the output). Solving Linear Systems: Iterative Methods For very large systems—particularly sparse systems where the matrix contains mostly zeros—direct methods become impractical. Iterative methods generate a sequence of approximations that gradually converge to the solution. Jacobi iteration and Gauss-Seidel iteration are straightforward methods that repeatedly refine each component of the solution using current estimates of the other components. Successive over-relaxation (SOR) accelerates convergence by extrapolating the improvements found in Gauss-Seidel. Conjugate gradient method is particularly effective for symmetric positive-definite matrices and often converges faster than the previous methods. These methods are memory-efficient (you don't store the full matrix, only how to compute $Ax$ for any vector $x$) and naturally parallelize for computation on multiple processors. Solving Nonlinear Equations When your equation isn't linear—for example, $x^3 - 2x + 1 = 0$—you need root-finding algorithms. Newton's method is popular when you know the derivative of $f(x)$. It starts with an initial guess $x0$ and iteratively improves it using: $$x{n+1} = xn - \frac{f(xn)}{f'(xn)}$$ The method converges rapidly (quadratic convergence) when started near the true root, but can diverge if started far away or if the derivative is small. Other methods like the bisection method are slower but more robust for difficult cases. Linearization Linearization converts a nonlinear problem into a linear one by approximating the nonlinear function with its linear (first-order) Taylor approximation around a known point. While this introduces error, it allows you to use the fast, well-understood methods for linear systems. This approach works well when the nonlinearities are weak or when you only need an approximate solution. Eigenvalue and Singular Value Problems Some problems are fundamentally about understanding the structure and properties of matrices themselves. Eigenvalue Decomposition An eigenvalue $\lambda$ and corresponding eigenvector $v$ of a matrix $A$ satisfy $Av = \lambda v$. Geometrically, multiplying by $A$ just scales the eigenvector by the eigenvalue—no rotation occurs. Eigenvalue decomposition reveals the natural "directions" and scaling factors inherent in a matrix. <extrainfo> A practical application is spectral image compression. Large images can be represented as matrices. Eigenvalue decomposition identifies the most important components of the image (those with largest eigenvalues) and discards less important details, achieving compression while preserving the recognizable image structure. </extrainfo> Singular Value Decomposition The singular value decomposition (SVD) is a more general factorization: $A = U\Sigma V^T$, where $U$ and $V$ are orthogonal matrices and $\Sigma$ is a diagonal matrix of non-negative values called singular values. Unlike eigenvalue decomposition, SVD works on any matrix (not just square ones) and is more numerically stable. <extrainfo> SVD underlies principal component analysis (PCA), a statistical technique that identifies the directions of maximum variance in high-dimensional data. By keeping only the components corresponding to the largest singular values, PCA reduces dimensionality while preserving most of the important structure in the data. This is fundamental to machine learning and data analysis. </extrainfo> Optimization Problems Beyond simply evaluating or solving, optimization finds the "best" solution according to some criterion. Linear Programming Linear programming problems have linear objective functions and linear constraints. For example: "Maximize profit = 3x + 2y subject to x + y ≤ 10, x ≥ 0, y ≥ 0." The feasible region (the set of points satisfying all constraints) forms a polygon, and the optimal solution occurs at a vertex of this polygon. The simplex method exploits this structure. It starts at one vertex and moves to an adjacent vertex with a better objective value, continuing until no improvement is possible. Despite being exponential in worst-case complexity, it's remarkably efficient in practice and was revolutionary when developed. Constrained Optimization: Lagrange Multipliers Many real problems have constraints: maximize revenue subject to production limits, minimize cost subject to quality standards, etc. The method of Lagrange multipliers elegantly converts a constrained optimization problem into an unconstrained one. For a problem "optimize $f(x)$ subject to $g(x) = 0$," you instead solve the unconstrained problem of finding the critical points of the Lagrangian function: $$L(x, \lambda) = f(x) - \lambda g(x)$$ The multiplier $\lambda$ (called the Lagrange multiplier) has a useful interpretation: it tells you how much the optimal value would improve if the constraint were relaxed slightly. This makes Lagrange multipliers valuable even beyond just solving the optimization problem. Numerical Integration (Quadrature) Computing definite integrals analytically is often impossible, especially for functions known only at discrete points. Numerical integration (also called quadrature) approximates $\inta^b f(x)\,dx$ using function values at selected points. Newton-Cotes Formulas The simplest approach: divide the interval $[a,b]$ into small subintervals and use simple geometric shapes to approximate the area under the curve. The midpoint rule approximates each subinterval with a rectangle of height equal to the function value at the interval's midpoint. Simpson's rule fits a parabola through three consecutive points and integrates the parabola exactly. It's remarkably accurate despite its simplicity: it exactly integrates all cubic polynomials (not just quadratics), giving it higher accuracy than expected. Gaussian Quadrature Rather than using equally-spaced evaluation points (as Newton-Cotes formulas do), Gaussian quadrature optimizes the choice of both the evaluation points and the weights assigned to each function value. This optimization dramatically improves accuracy. For a given number of function evaluations, Gaussian quadrature achieves higher accuracy than Newton-Cotes formulas. The trade-off: the optimal points depend on the interval and function class (they're specially tabulated), whereas Newton-Cotes always uses evenly-spaced points. High-Dimensional Integration In high dimensions, traditional quadrature becomes impractical—the number of points needed grows exponentially with dimension. For these problems: Monte Carlo integration randomly samples points and computes the average function value, scaled by the region's volume. Simple but converges slowly. Quasi-Monte Carlo integration uses carefully-chosen non-random sequences of points that spread more evenly than random samples, improving convergence rates. Sparse-grid methods use clever combinations of one-dimensional quadrature rules to handle modest dimensions more efficiently than full tensor products. Differential Equations Many physical phenomena are modeled by differential equations describing rates of change. Numerical methods solve these when analytical solutions are unavailable. Ordinary Differential Equations An ordinary differential equation (ODE) involves a function of one variable and its derivatives. The simplest is the initial value problem: given $\frac{dy}{dt} = f(t, y)$ and an initial condition $y(t0) = y0$, find $y(t)$ for $t > t0$. The Euler method is the simplest numerical solver. It advances the solution step-by-step: $$y{n+1} = yn + h f(tn, yn)$$ where $h$ is the step size and the term $h f(tn, yn)$ estimates how much $y$ changes over time interval $h$ using the estimated slope. Euler's method is intuitive and easy to implement, but only moderately accurate. Better methods (Runge-Kutta methods) use evaluations at multiple points within each step to achieve higher accuracy. Partial Differential Equations A partial differential equation (PDE) involves functions of multiple variables (like space and time) and their partial derivatives. Solving PDEs analytically is often impossible. Numerical solutions discretize the continuous domain into a grid and approximate derivatives using differences between nearby grid points. Three main approaches: Finite difference method replaces derivatives with difference quotients: $\frac{\partial u}{\partial x} \approx \frac{u(x+h) - u(x)}{h}$. Finite element method divides the domain into small elements, approximates the solution as a simple function (like a polynomial) on each element, and enforces continuity at element boundaries. Particularly useful for irregular domains. Finite volume method divides the domain into control volumes and ensures conservation laws (like conservation of mass) hold within each volume. Naturally suited to conservation laws. All three approaches convert the continuous PDE into a system of algebraic equations at the grid points, which you then solve using the linear system methods discussed earlier.