Numerical analysis

Numerical analysis studies algorithms for approximating mathematical problems and controlling their errors.

A finite volume method evolves cell integrals or cell averages by differences of shared interface numerical fluxes. For a scalar conservation law, $U_j^{n+1}=U_j^n-(k/d)(F_{j+1/2}-F_{j-1/2})$. Sharing the same flux between neighboring cells makes conservation follow by telescoping.

A scalar finite-volume update is monotone if increasing any input state cannot decrease any output state. A consistent conservative update preserves constants and has an invariant interval by comparison. On a periodic grid, or an infinite grid with summable perturbations, conservation and monotonicity imply L1 contraction of a monotone conservative scheme. Translation invariance then implies the total variation diminishing scheme property.

A translation-invariant conservative monotone update decreases discrete total variation. If $S$ is the one-cell shift, $TS=ST$, and L1 contraction of a monotone conservative scheme applied to $U,SU$ gives ${\sum }_j|(TU{)}_{j+1}-(TU{)}_j|\le {\sum }_j|U_{j+1}-U_j|$. This prevents creation of new total variation, though it does not make every discontinuity sharply resolved.

Let $T$ be a monotone conservative scalar update on a periodic grid. Since $T(U\vee V)\ge TU,TV$, one has $(TU-TV{)}_{+}\le T(U\vee V)-TV$. Conservation makes the sum of the right side equal to $\sum (U-V{)}_{+}$. Reverse $U,V$ and add to obtain $\Vert TU-TV{\Vert }_1\le \Vert U-V{\Vert }_1$. On an infinite grid the same proof holds for summable differences with a Lipschitz flux on the state range.

The Engquist-Osher method is the finite volume method using the Engquist-Osher flux. Under $(k/d)\sup |f\prime |\le 1$ on the invariant state interval, its update is monotone. It preserves the state range, is contractive for summable perturbations in $L^1$, and decreases discrete total variation.

The Engquist-Osher numerical flux is $F(a,b)=f(0)+{\int }_0^a\max (f\prime (s),0)ds+{\int }_0^b\min (f\prime (s),0)ds$. It sends positive-speed contributions from the left state and negative-speed contributions from the right state. Its consistency $F(u,u)=f(u)$ and opposite monotonicities in the two arguments yield a monotone conservative scheme under the appropriate Courant–Friedrichs–Lewy condition. For the Inviscid Burgers equation, it is $\frac{1}{2}\max (a,0{)}^2+\frac{1}{2}\min (b,0{)}^2$.

A numerical flux approximates the transported quantity through a cell interface using nearby states. For a two-state flux $F(a,b)$ approximating a scalar conservation law with flux $f$, consistency requires $F(u,u)=f(u)$. A monotone conservative scheme uses a flux nondecreasing in its left state and nonincreasing in its right state, with a time step making the complete update monotone.

For a continuous increasing scalar function, retain an interval whose endpoints bracket a target value, evaluate its midpoint, and keep the half that still brackets the target. Each step halves the interval length. Applied to a strictly increasing water-level equation, this finds its unique solution to any specified error tolerance.

An implicit time-stepping method determines a new solution value through an equation involving that same value. A small-step contraction mapping argument often establishes local solvability when the vector field is Lipschitz continuous. Examples include the Backward Euler method, implicit Runge-Kutta methods and implicit linear multistep methods. Solvability of the update is a separate issue from stability of a numerical method.

A multigrid method combines relaxation of high-frequency fine-grid errors with coarse-grid correction of smooth errors. Residual restriction, an approximate coarse error solve, error prolongation and post-smoothing form a two-grid cycle; recursive coarse solves give multilevel cycles.

For residual $r_h=b_h-A_h{u}_h$, restrict to a coarse grid, solve $A_H{e}_H=R{r}_h$ approximately, and update $u_h$ by $P{e}_H$. This corrects error that is smooth on the fine grid but representable on the coarse grid.

For an initial value problem, convergence means that the largest global error on each fixed finite time interval tends to zero as the step size tends to zero and the starting errors vanish. For linear multistep methods, the Dahlquist equivalence theorem separates this requirement into consistency of a numerical method and zero-stability.

The global error is the difference between the computed and exact solutions after accumulated numerical steps, rather than the local truncation error of a single step started from exact data. An estimate for the stability of a numerical method controls how initial errors and local defects accumulate. For a stable order-$p$ time integrator with starting errors $O(h^p)$, local defects $O(h^{p+1})$ accumulated over $O(1/h)$ steps give global error $O(h^p)$ on a fixed time interval.

The midpoint quadrature rule approximates an integral on $[a,b]$ byIt is exact for linear polynomials.

A truncation error is the error introduced by replacing an exact mathematical operation with a finite approximation, such as replacing a derivative by a finite difference.

A spectral method approximates a differential equation in a global basis and enforces the equation after projection or collocation. Smooth and analytic solutions often give rapidly convergent coefficients.

A Fourier spectral method represents a periodic unknown by finitely many Fourier modes and converts differentiation into multiplication by wavenumber.

A spectral truncation retains a finite set of basis modes and sets all omitted coefficients to zero.

For the two-periodic problemthe Fourier coefficients satisfy the tridiagonal system

The Fourier–Galerkin method seeks a finite Fourier series whose residual is orthogonal to every retained Fourier mode.

For a numerical evolution, stability means that perturbations of the initial data and forcing remain bounded over a fixed time interval, with a bound independent of the discretization parameters. For a linear homogeneous recurrence, a mesh-uniform bound on powers of its update matrix provides this estimate. Von Neumann stability analysis and the energy method are two complementary ways to establish it.

A numerical method is consistent when its local defect tends to zero at the rate required to approximate the differential equation. Stability then determines whether these local defects accumulate into a convergent global approximation.

For a time-stepping formula written with an unscaled step residual $R_h$, order at least $p$ means that substituting a smooth exact solution gives $R_h=O(h^{p+1})$. Under the relevant stability of a numerical method assumptions and starting errors $O(h^p)$, this produces global error $O(h^p)$. Residuals divided by the step size instead use the equivalent convention $O(h^p)$.

The energy method differentiates a norm or quadratic energy along a differential equation or numerical evolution and estimates the resulting terms. Dissipative, skew-adjoint, and boundary contributions can yield stability or conservation without diagonalizing the operator.

For the centered second difference and centered first difference on a finite interval with zero Dirichlet boundary conditions, the diffusion matrix is symmetric negative definite and the drift matrix is skew-symmetric. Discrete summation by parts gives the displayed mesh-weighted identity. It proves contraction in the discrete L2 norm without needing periodic modes, positive stencil coefficients or simultaneous diagonalization.

For a zero-boundary centered semidiscretization $U^{\prime }=D_{h}U+V_{h}U$, where $D_h$ is the scaled Dirichlet discrete Laplacian and $V_h=\mathrm{diag}(a_j)$ with $a_j\le a_{+}$, summation by parts gives $(U,D_{h}U{)}_h\le 0$. Hence $\frac{1}{2}(d/dt)\Vert U{\Vert }_h^2\le a_{+}\Vert U{\Vert }_h^2$. The Gronwall inequality yields $\Vert U(t){\Vert }_h\le e^{a_{+}t}\Vert U(0){\Vert }_h$, a mesh-independent stability bound on every fixed finite interval. Positive reaction coefficients can permit physical growth without violating this form of stability.

For the weak energy solution of a variable-coefficient wave equation, uniqueness can be proved without using $u_t$ as a spatial $H_0^1$ test function. For a difference with zero data, fix $s$ and test by $v(t)={\int }_t^{s}u(r)dr$ for $t<s$, extended by zero afterwards. Put $z(t)={\int }_0^{t}u(r)dr$, so $v(t)=z(s)-z(t)$. Integration by parts in time for the symmetric principal form, and in space for the first-order terms, giveFor short enough intervals the last term is absorbed. The Gronwall inequality forces $u=0$, and repetition proves uniqueness on the full interval. This antiderivative test is valid at the stated space-time $H^1$ regularity, unlike an unqualified direct energy test by $u_t$.

The finite element method approximates a variational problem in a finite-dimensional space of piecewise polynomial functions and determines the coefficients by testing against the same trial space.

The spatial Galerkin method for the heat equation gives a mass matrix $M$ and a diffusion stiffness matrix $S$. With zero forcing, $\frac{d}{dt}(c^{T}Mc)=-2c^{T}Sc\le 0$. The quadratic form $c^{T}Mc$ is the squared L2 norm of the represented finite-element function, proving an energy contraction. For nonzero forcing, the corresponding identity has the additional term $2c^{T}F$, which can be estimated using the Cauchy-Schwarz inequality.

A finite element specifies a reference region, a finite-dimensional space of local trial functions, and degrees of freedom that uniquely determine each such function. Assembling compatible copies over a finite element mesh gives a global trial space. For a one-dimensional linear element, the region is an interval, the functions are affine functions, and the degrees of freedom are the endpoint values.

A conforming finite element space is a finite-dimensional vector subspace of the solution's variational space $V$. For a second-order problem with homogeneous Dirichlet boundary conditions, continuous piecewise-polynomial functions vanishing on the boundary form such a subspace of the zero-boundary Sobolev space $H_0^1(\Omega )$. Conformity allows exact Galerkin orthogonality and the Céa lemma.

A finite element mesh partitions the domain into simple elements, such as intervals, triangles or tetrahedra. Piecewise-polynomial trial functions and local integration assemble the stiffness matrix and mass matrix.

A family of simplicial meshes is shape-regular when the ratio of each element's diameter to its inscribed-ball radius stays uniformly bounded. This excludes degenerating shapes and makes finite element interpolation estimates uniform in the mesh size. The element sizes may still vary substantially across the domain.

On a shape-regular mesh in dimensions at most three, nodal piecewise-linear interpolation of $u\in H^2(\Omega )$ satisfieswhere $h$ is the largest element diameter and $C$ is independent of $h$. For homogeneous Dirichlet boundary conditions, $I_{h}u$ belongs to the conforming finite element space. Scaling each element to a reference element gives the estimate; shape regularity bounds the scaling constants, while subtraction of an affine function leaves an error controlled by the second derivatives.

The mass matrix is the Gram matrix of a finite element basis in the L2 inner product. Its coefficient quadratic form is the squared L2 norm of the represented trial function, so it is a positive-definite matrix for a linearly independent basis.

A mass eigenstate diagonalizes the Hermitian mass operator or mass matrix and has a definite mass. In a flavor basis with off-diagonal mass terms, it is a linear combination of flavor states. For unstable particles, decay propagation instead involves an effective complex matrix; its eigenvectors need not inherit the orthogonality of a Hermitian mass problem.

On a one-dimensional uniform finite element mesh of spacing $h$, let $w$ be an affine function and ${\varphi }_i$ the interior piecewise-linear hat functions. The weighted mass matrix has entriesand zero entries for nonadjacent indices. The diagonal product ${\varphi }_i^2$ is symmetric about $x_i$ and integrates to $2h/3$; the neighbouring product is symmetric about $(x_i+x_{i+1})/2$ and integrates to $h/6$. The odd part of an affine function integrates to zero around each centre, giving the formulas. This evaluates a variable reaction term exactly without substituting a constant coefficient or a quadrature approximation.

The stiffness matrix represents a variational bilinear form $a$ in a finite element basis. A symmetric coercive bilinear form gives a positive-definite matrix. Local support of the piecewise-linear hat functions produces a sparse matrix, allowing efficient numerical solution.

For $-u^{\prime \prime }+u=f$ with homogeneous Dirichlet boundary conditions, the conforming piecewise-linear hat functions on a uniform grid give a diffusion stiffness matrix $S$ and a mass matrix $M$. Their nonzero entries are $S_{ii}=2/h$, $S_{i,i+1}=-1/h$, $M_{ii}=2h/3$, and $M_{i,i+1}=h/6$, with symmetric counterparts. The sum is a positive-definite matrix because its coefficient quadratic form is $\int ((u_h^{\prime }{)}^2+u_h^2)$. For forcing $f(x)=x$, the load at $x_i$ is $h{x}_i$.

For periodic $u_t=u_x$, the conforming Galerkin method has $M\dot{U}=CU$, with mass matrix $M_{ij}=\int {\varphi }_i{\varphi }_j$ and $C_{ij}=\int {\varphi }_i{\varphi }_j^{\prime }$. Integration by parts and periodicity give $C^T=-C$. ThereforeThe conserved quantity is the squared L2 norm of the finite element function. On a uniform one-dimensional mesh of spacing $h$, the piecewise-linear hat functions give $M_{ii}=2h/3$, $M_{i,i\pm 1}=h/6$, $C_{i,i+1}=1/2$ and $C_{i,i-1}=-1/2$.

A cubic Hermite finite element is piecewise cubic and joins with continuous first derivative. It is therefore conforming for fourth-order variational problems whose energy space lies in $H^2$.

The Ritz method minimizes a quadratic energy over a finite-dimensional trial space. For basis functions ${\varphi }_j$ and bilinear form $a$, its linear system is $Ka=F$ with $K_{jk}=a({\varphi }_k,{\varphi }_j)$ and $F_j=⟨f,{\varphi }_j⟩$.

For a bounded symmetric coercive bilinear form on a Hilbert space, a solution of $a(u,w)=\ell (w)$ minimizes $J(v)=a(v,v)/2-\ell (v)$. The displayed identity proves strict minimality and uniqueness. Restricting to a conforming finite element space makes the discrete Ritz method equivalent to the Galerkin method. Galerkin orthogonality gives exact best approximation in the energy norm; the Céa lemma translates it into a reference-norm estimate.

A symmetric coercive bilinear form $a$ defines the energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$. In this norm the Ritz method is an orthogonal projection onto its finite-dimensional trial space.

If $u$ solves a variational problem and $u_h$ is its conforming Galerkin method approximation in $V_h$, thenfor every $v_h\in V_h$.

For a continuous coercive bilinear form with continuity constant $M$ and coercivity constant $\alpha$, a conforming finite-element solution satisfiesFor a symmetric problem measured in its energy norm, the constant is one.

For any $v_h\in V_h$, Galerkin orthogonality gives $a(e,e)=a(e,u-v_h)$, where $e=u-u_h$. Coercivity and continuity yield $\alpha \Vert e{\Vert }_V^2\le M\Vert e{\Vert }_V\Vert u-v_h{\Vert }_V$. Divide by $\Vert e{\Vert }_V$ when it is nonzero and take the infimum. In the symmetric case, the Cauchy-Schwarz inequality in the energy norm gives the constant one.

For a symmetric elliptic variational problem with Galerkin orthogonality, let $e=u-u_h$ and solve the dual problem $a(v,z)=(e,v{)}_{L^2}$ for every $v\in H_0^1(\Omega )$. If elliptic regularity gives $\Vert z{\Vert }_{H^2}\le C\Vert e{\Vert }_{L^2}$, then a finite element interpolation estimate givesThus $\Vert e{\Vert }_{L^2}\le Ch\Vert e{\Vert }_{H^1}$. Combining this with the Céa lemma and an $O(h)$ energy norm error gives an $O(h^2)$ L2 norm error. The additional order depends on the dual regularity assumption, which can fail on unsuitable domains.

A piecewise-linear hat function is one at one mesh node, zero at all other nodes, and affine on each adjacent mesh interval. Interior hat functions form the standard nodal basis for continuous piecewise-linear finite elements with homogeneous Dirichlet boundary conditions.

An arithmetic operation is one scalar addition, subtraction, multiplication, or division. Counting such operations gives a basic model of an algorithm's computational cost.

For a real monic quadratic $z^2+a_{1}z+a_2$, both roots lie in the closed unit disk whenwith the usual simplicity requirement for roots on the unit circle.

For a complex polynomial $p(z)=a{z}^2+bz+c$, assume $|a|>|c|$ and put $D=|a{|}^2-|c{|}^2$ and $E=\overline{a}b-c\overline{b}$. Both roots of a polynomial lie in the closed unit disk if and only if $|E|\le D$. Strict inequality puts both roots in the open disk. Under the strict product assumption $|c/a|<1$, any unit-modulus root is automatically simple.

To prove the strict version, define $p^{\#}(z)=z^2\overline{p(1/\overline{z})}$ and $q=\overline{a}p-c{p}^{\#}=z(Dz+E)$. On the unit circle $|p^{\#}|=|p|$, so the Rouche theorem shows that a stable $p$ and $q$ have the same two interior zeros. Conversely $Dp=aq+c{q}^{\#}$, and the same Rouche theorem implies that two interior zeros of $q$ give two interior zeros of $p$. The nonzero root of $q$ is $-E/D$, proving the criterion. The closed version follows by continuity, replacing $b$ by $(1-\epsilon )b$ in the sufficiency direction. In the necessity direction with a boundary root, replace $p(z)$ by $p(z/r)$ for $r<1$ close to one and let $r\to 1$. The product of the two roots of a polynomial excludes two unit-modulus roots and excludes a repeated unit-modulus root.

For bounded matrices $A_j$, the Lie product formula gives operator norm convergenceFor two summands, expanding their matrix exponentials gives the one-step difference below. For Hermitian matrices multiplied by $-i$, the factors are unitary operators and the first-order unitary product-formula error bound gives an explicit error estimate.

For bounded matrices, differentiating the symmetrized product gives $F^{\prime }-(A+B)F=\frac{1}{2}([e^{tB},A]e^{tA}+[e^{tA},B]e^{tB})$. The variation-of-constants formula turns this into an integral expression for $F-e^{t(A+B)}$. This retains the noncommutativity of the factors explicitly; cancellation of their leading opposite commutators explains why the averaged product is more accurate than either ordering alone.

For bounded matrices,Commutativity removes this splitting error, though first-order substeps such as Backward Euler method retain their own quadratic local errors.

For Hermitian matrices $H_1,\dots ,H_m$ and $t\ge 0$, let $S(t)={\prod }_{j=1}^m{e}^{-it{H}_j}$. ThenHere every norm is the spectral norm. A short proof starts with two summands. Differentiate $F(s)=e^{-i(t-s)(A+B)}{e}^{-isA}{e}^{-isB}$; its derivative has norm at most $\Vert [e^{-isA},B]\Vert$. Differentiating $e^{-iuA}B{e}^{iuA}$ and integrating yields $\Vert [e^{-isA},B]\Vert \le s\Vert [A,B]\Vert$, because unitary operators preserve the norm. Integrating $s$ from zero to $t$ gives $t^2\Vert [A,B]\Vert /2$. Inductively separate $H_1$ from the remaining sum, use the triangle inequality on their commutator, and apply the telescoping bound for products of operators to obtain the displayed many-term bound.

Repeating steps of size $t/k$ and telescoping across $k$ steps givesFor a chain of $O(n)$ bounded nearest-neighbor terms, only $O(n)$ pairs fail to commute. Consequently first-order product-formula Hamiltonian simulation has error $O(n{t}^2/k)$ and uses $O(nk)$ constant-size gates. The general bound is also given in Proposition 9 of Childs and collaborators' analysis of Trotter error.

Strang splitting approximates $e^{h(A+B)}$ by $e^{hA/2}{e}^{hB}{e}^{hA/2}$. Its time symmetry cancels the quadratic local commutator defect, giving local error $O(h^3)$ and global order two.

If real matrices satisfy $A^T=A\le 0$ and $B^T=-B$, then $e^{tA}$ is a contraction and $e^{tB}$ is an orthogonal matrix for $t\ge 0$. Hence $\Vert e^{kA/2}{e}^{kB}{e}^{kA/2}{\Vert }_2\le 1$ for every $k\ge 0$. The bound is independent of commutativity and gives unconditional stability for a centered convection-diffusion semidiscretization with homogeneous Dirichlet boundary conditions. It does not itself bound the mesh-dependent commutators that enter a time-error estimate.

If $S_h$ is a symmetric second-order splitting, then $S_{\gamma h}{S}_{\delta h}{S}_{\gamma h}$ is fourth order for $\gamma =(2-2^{1/3}{)}^{-1}$ and $\delta =-2^{1/3}(2-2^{1/3}{)}^{-1}$. The identities $2\gamma +\delta =1$ and $2{\gamma }^3+{\delta }^3=0$ preserve the total step and cancel the leading odd defect.

For $f\in C^2[a,b]$, the composite midpoint rule with mesh width $h$ has local error $O(h^3)$ and global error $O(h^2)$.

If $f(x)=x^{-1/2}g(x)$ with $g$ smooth near zero, the substitution $x=s^2$ givesThe transformed integrand is smooth, so an ordinary second-order midpoint rule recovers its usual convergence rate.

Given distinct nodes $x_0,\dots ,x_n$ and values $f(x_i)$, there is a unique polynomial of degree at most $n$ taking those values.

An interpolation node is one of the distinct argument values at which an interpolating function is required to match prescribed data.

The interpolating polynomial is

For distinct interpolation nodes $x_0,\dots ,x_n$, the polynomial ${\ell }_i(x)={\prod }_{j\ne i}(x-x_j)/(x_i-x_j)$ has ${\ell }_i(x_j)={\delta }_{ij}$. These cardinal polynomials form a basis of the polynomials of degree at most $n$, and give $p(x)={\sum }_{i}p(x_i){\ell }_i(x)$.