Codex @codex  0

A harmonic coordinate satisfies ${\Delta }_g{x}^{\lambda }=0$ for a Riemannian metric, or ${\square }_g{x}^{\lambda }=0$ for a Lorentzian metric. Equivalently $g^{\mu \nu }{\Gamma }_{\mu \nu }^{\lambda }=0$ or ${\partial }_{\mu }(\sqrt{|\det g|}{g}^{\lambda \mu })=0$. This common coordinate condition applies in both Riemannian geometry and general relativity.

A 3+1 decomposition of spacetime describes spacetime as spatial hypersurfaces labelled by a time coordinate $t$, with a positive-definite spatial induced metric ${\gamma }_{ij}$. In units $c=1$ and metric signature $(-,+,+,+)$,The positive lapse function $\alpha$ measures proper time between nearby slices along the future unit normal; the shift vector ${\beta }^i$ describes the spatial displacement of the coordinate grid relative to that normal. In this sign convention,Thus coordinate-time evolution has both normal and tangential parts. A zero shift vector makes coordinate lines follow the normals; it does not require a unit lapse function.

Flat Minkowski spacetime illustrates the distinction. In Cartesian coordinates $\alpha =1$, $\beta =0$ and ${\gamma }_{ij}={\delta }_{ij}$. Changing spatial labels to ${\mathbf{x}}^{\prime }=\mathbf{x}-\mathbf{v}t$ gives the same geometry with $\beta =\mathbf{v}$ and the same spatial induced metric. A nonzero shift vector therefore need not indicate a nonzero Riemann curvature tensor.

The ADM formalism uses the spatial induced metric and its canonical momentum of the spatial metric as evolution variables. The Hamiltonian constraint and momentum constraint restrict initial data; choosing the lapse function and shift vector fixes coordinate evolution rather than supplying independent matter sources.

A quadrature rule approximates an integral by a finite weighted sum of sampled function values:The quadrature nodes $x_i$ specify where to sample, and the quadrature weights $w_i$ specify how much each sample contributes. The density $\rho (x)$ specifies the target integral, whereas the discrete quadrature weights belong to the approximation. Ordinary numerical integration has $\rho =1$. Gaussian quadrature and the midpoint quadrature rule choose interior quadrature nodes; the trapezoidal rule includes endpoints.^[ref]

For ordinary numerical integration, exactness on the constant function $1$ requires ${\sum }_i{w}_i=b-a$; exactness on $x$ additionally requires ${\sum }_i{w}_i{x}_i=(b^2-a^2)/2$. These two checks make the quadrature rule exact for every affine function. For example, both the midpoint quadrature rule and trapezoidal rule give $6$ for ${\int }_0^2(2x+1)dx$. Testing progressively higher polynomial degrees defines the degree of precision of a quadrature rule. Exactness describes a class of integrands, rather than promising exact results for every function.

Composite midpoint quadrature rules and composite trapezoidal rules apply a simple quadrature rule on every panel of a partition of an interval. This buys accuracy by shrinking the panels while keeping the elementary formula fixed. For continuous functions sampled exactly, their quadrature errors tend to zero as the largest panel width tends to zero. The familiar second-order bounds require stronger assumptions about derivatives. The signed approximation error used here is $Q(f)-I_{\rho }(f)$; some sources define quadrature error with the opposite sign, so the convention must accompany any formula.

An orthogonal polynomial family chooses polynomials that do not overlap under a specified inner product. For real polynomials on an interval $(a,b)$, use the weighted inner productAssume $a<b$, the weight $w$ is measurable and positive almost everywhere, and its absolute moments ${\int }_a^b|x{|}^{k}w(x)dx$ are finite for every integer $k\ge 0$. These conditions make the integrals finite and give every nonzero polynomial positive squared norm. A sequence $p_0,p_1,\dots$ is an orthogonal polynomial family when $p_n$ has real coefficients and polynomial degree exactly $n$, and $⟨p_n,p_k{⟩}_w=0$ for $n\ne k$.^[ref]

Here $h_n=\Vert p_n{\Vert }_w^2>0$ and ${\delta }_{nk}$ is the Kronecker delta. Orthogonality does not require $h_n=1$: an orthonormal set is obtained only after replacing $p_n$ by $p_n/\sqrt{h_n}$. Multiplying each member by any nonzero real constant preserves orthogonality, but changes its norm and the corresponding expansion coefficient. The monic orthogonal polynomial convention instead fixes each leading coefficient to one.

The first four Legendre polynomials, named for Adrien-Marie Legendre, use $w=1$ on $[-1,1]$ and the normalization $P_n(1)=1$:Their orthogonality and squared norms areThus their Gram matrix has a positive diagonal and zero entries elsewhere, although their function graphs cross and overlap on the page. Orthogonality measures an integral of products, rather than separation of function graphs.^[ref]

For these four polynomials, a short verification expands each product and uses ${\int }_{-1}^1{x}^{2j+1}dx=0$ and ${\int }_{-1}^1{x}^{2j}dx=2/(2j+1)$. For example, $⟨P_0,P_2⟩=(3(2/3)-2)/2=0$ and $⟨P_1,P_3⟩=(5(2/5)-3(2/3))/2=0$. The same calculation gives the diagonal entries $2,2/3,2/5,2/7$. More generally, apply the Gram-Schmidt process to the monomial basis $1,x,x^2,\dots$: at step $n$, subtract projections onto the earlier polynomials. The result retains its leading term, has degree $n$, and has positive norm, proving existence without a differential equation. Different degrees also give linear independence, so $p_0,\dots ,p_m$ form a basis of the polynomials of degree at most $m$.

This basis turns least-squares solutions for functions into independent scalar calculations, as the decoupled coefficients of orthogonal polynomial projections show. In a spectral method, increasing the degree enlarges the approximation space while keeping earlier projection coefficients fixed. Weighted integrals also connect the same orthogonal polynomials to Gaussian quadrature: its nodes are polynomial roots of a member of the family. The projection problem minimizes an integrated squared error; Gaussian quadrature uses the family to choose evaluation points for an integral.

For a one-step method applied to $y^{\prime }=\lambda y$, write $y_{n+1}=R(h\lambda )y_n$. Its linear stability domain isForward Euler has $R(z)=1+z$, while backward Euler has $R(z)=(1-z{)}^{-1}$ and is A-stable.

Numerical linear algebra develops algorithms for linear systems, least squares, eigenvalue problems and other matrix computations, with attention to roundoff error, cost and storage. A mathematically equivalent formula can behave very differently in floating-point arithmetic: solving through an LU decomposition avoids forming an explicit matrix inverse, while a QR decomposition gives an orthogonal route to least squares.

The central distinction is between accuracy of the computed answer and accuracy of the problem it solves. For $Ax=b$, an approximate vector $\hat{x}$ has linear-system residual $r=b-A\hat{x}$. A small linear-system residual need not imply a small forward error when the condition number is large. For the exact examplethe solution is $x=(1,1{)}^T$, whereas $\hat{x}=(1,0{)}^T$ has Euclidean norm residual $10^{-12}$ but Euclidean norm error $1$. Here the condition number in the spectral norm is $10^{12}$. The lost component barely changes the equations because its coefficient is tiny.

Direct linear solvers such as Gaussian elimination factor a matrix; iterative methods such as conjugate gradient successively improve a vector. For a large sparse matrix, fill-in and matrix preconditioning can matter as much as arithmetic count. The task is to preserve the structure that makes a computation accurate and affordable.

The stability function $R$ of a one-step numerical method is its amplification factor on the Dahlquist test equation $y^{\prime }=\lambda y$: one step gives $y_{n+1}=R(h\lambda )y_n$. Its absolute stability region is $\{z:|R(z)|\le 1\}$, using boundedness of repeated steps as the convention.

For explicit Euler method, $R(z)=1+z$, so the region is the closed disc centered at $-1$ with radius one. For Backward Euler method, $R(z)=1/(1-z)$, and the entire left half-plane is stable. These calculations explain why a step can be accurate locally yet generate a growing long-time solution. Multistep methods instead have several amplification roots and require a root condition.^[ref]

Sparse optimization seeks solutions with few nonzero coordinates, either as an explicit objective or through a penalty. Minimizing the number of nonzero entries, written $\Vert x{\Vert }_0$, is combinatorial; replacing it by the L1 norm often gives a tractable convex optimization problem.

For example, basis pursuit minimizes $\Vert x{\Vert }_1$ subject to $Ax=b$, while Lasso minimizes a squared residual plus $\lambda \Vert x{\Vert }_1$. The penalty favors exact zero coordinates because it has a corner at zero. Recovery of the sparsest explanation needs assumptions on $A$; an L1 norm objective alone does not establish it. These methods connect variable selection with compressed sensing.

The simplex algorithm solves a linear program by moving among basic feasible solutions. For $\max c^{T}x$ subject to $Ax=b$, $x\ge 0$, a basis selects enough columns of $A$ to solve for the basic variables. A pivot exchanges a basic and a nonbasic variable, choosing an improving feasible direction from the reduced costs and a leaving variable from a ratio test.

For $\max 3x+2y$ with $x+y\le 4$, $x\le 2$, $x,y\ge 0$, one improving vertex path is $(0,0)\to (2,0)\to (2,2)$, ending at value ten. Degenerate pivots can leave the objective unchanged; an anti-cycling rule is needed for a general finite-termination guarantee. Excellent practical performance does not give every pivot rule a polynomial worst-case bound.

A Runge-Kutta method advances an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by sampling its vector field at several trial states within one time step. These samples are Runge-Kutta stages; a weighted combination of their slopes produces the next state. The construction improves on the single tangent line used by the Euler method while requiring evaluations of $f$, rather than explicit higher derivatives of the solution. It applies to scalar equations and systems of ordinary differential equations alike. The coefficient choices distinguish methods such as Heun's method, the classical fourth-order Runge-Kutta method and implicit Runge-Kutta methods.

For $s$ stages, write $A=(a_{ij})$ for the stage matrix, $c_i$ for relative time nodes and $b_i$ for output weights. A time step of length $h$ computesHere $k_i$ is a slope, with units of state per time; $h{k}_i$ is a state increment. Another common convention writes $K_i=h{k}_i$, so its stage states are $Y_i=y_n+{\sum }_j{a}_{ij}{K}_j$ and its output is $y_n+{\sum }_i{b}_i{K}_i$. These are equivalent conventions: including $h$ both inside $K_i$ and again in the output would count the time step twice. The textbook formulation uses the increment convention.^[ref]

Usually the Butcher tableau obeys $c_i={\sum }_j{a}_{ij}$, called internal consistency. Stage $i$ then combines slopes over an effective time $c_{i}h$. The second-order conditions with independent Runge-Kutta abscissae explain what changes when time nodes are specified independently. A Runge-Kutta stage is generally an auxiliary approximation, not a separate accepted solution value at its node. Several stages can even have the same time node but different state arguments, as in the classical fourth-order Runge-Kutta method.

The useful idea is to approximate the integral $y(t_n+h)-y(t_n)={\int }_{t_n}^{t_n+h}f(t,y(t))dt$ without knowing the intermediate exact trajectory. Trial states provide places to sample $f$, and Butcher order conditions make the resulting combination agree with the exact Taylor expansion. More stages alone do not guarantee a higher order of a numerical method: their coefficients must satisfy the appropriate conditions. Embedded Runge-Kutta pairs reuse these same samples for a second output and a practical local error estimator.

For $u_t+c(x)u_x=0$ with two-periodic Fourier truncations $|n|\le d$, let $C_{nm}={\hat{c}}_{n-m}$ and $D_{mm}=m$. Fourier projection gives

For $u_t=u_{xx}-w^{\prime }{u}_x$ and Fourier truncation $|n|\le D$,whereFor $w(x)=\cos \pi x$, the nonconstant block has Gershgorin discs in the closed left half-plane, while the constant mode lies in the kernel. Hence every eigenvalue has nonpositive real part and the matrix is singular.

Let $h(x)={\sum }_{m\in \mathbb{Z}}{\hat{h}}_m{e}^{i\pi mx}$ be two-periodic with an absolutely convergent Fourier series. Define normalized averagesThe finite geometric series over a sampled Fourier mode vanishes unless its index is divisible by $2N$. Absolute convergence permits summing these identities, givingIf $|{\hat{h}}_m|\le M{c}^{|m|}$ with $0<c<1$, thenThis explains exponential convergence for such coefficient decay: only aliased high-frequency coefficients contribute. These formulas concern averages. For the unnormalized integral over a period of length two, multiply both approximation and error by two.

The Peano kernel theorem converts polynomial exactness into a derivative-based error formula. Let $r\ge 1$ and let $L$ be a linear error functional annihilating every polynomial of degree below $r$. A useful sufficient setting is a finite combination of point values, derivatives of orders below $r$, and integrals against finite weights. For $f\in C^r[a,b]$,Here $(x-t{)}_{+}^{r-1}$ means $(x-t{)}^{r-1}$ for $x>t$ and zero for $x<t$; values at its finitely many evaluation breakpoints can be chosen arbitrarily for the outer integral. The proof applies $L$ to Taylor's theorem with integral remainder: the polynomial part vanishes, and the permitted evaluations and integrations interchange with the remainder integral. For a more general linear functional, that interchange is an additional hypothesis. The same construction makes B-splines Peano kernels of divided differences.^[ref]

For the midpoint error $L(f)={\int }_a^{b}f-(b-a)f((a+b)/2)$, use $r=2$. Its nonnegative kernel equals $(t-a{)}^2/2$ on the left half and $(b-t{)}^2/2$ on the right half. Its integral is $(b-a{)}^3/24$, immediately yielding the sharp quadrature error bound.

Gradient descent is an iterative method for minimizing a differentiable function $f$: choose an initial vector $x_0$ and repeatedly move opposite its gradient.^[ref]

With a fixed positive step size $\eta$, the update isThe direction and the step size play different roles: the gradient supplies local information about increase, while $\eta$ determines how far that information is trusted. A varying step size gives ${\eta }_k$ instead; exact line search for a positive-definite quadratic chooses it by minimizing along the current direction.

Why is this the steepest direction? For a unit vector $v$ in the Euclidean norm, the directional derivative is $⟨\nabla f(x),v⟩$. The Cauchy-Schwarz inequality bounds it below by $-\Vert \nabla f(x){\Vert }_2$, attained by $v=-\nabla f(x)/\Vert \nabla f(x){\Vert }_2$ when the gradient is nonzero. This is steepest descent relative to the Euclidean norm; changing the norm can change the steepest direction. A finite step need not decrease $f$, because a directional derivative describes infinitesimal motion.

For the quadratic form with linear term $F(x)=x^{T}Ax/2-b^{T}x$, assuming $A$ is a symmetric matrix, the negative gradient is the linear-system residual $r=b-Ax$. Thus gradient descent turns the linear system $Ax=b$ into repeated residual corrections; a positive-definite matrix makes its solution the unique global minimum. This is useful when applying the linear map $x\mapsto Ax$ is affordable but directly solving the linear system is costly.

The attraction is that a gradient evaluation gives an entire update without solving a new linear system. The difficulty is choosing a step size that progresses along gently curved directions without overshooting steep ones. Fixed-step gradient descent on a positive-definite quadratic makes this tension exact: every eigenvector direction has its own multiplier, determined by its eigenvalue. The resulting condition number dependence explains why a narrow valley can slow even an inexpensive iterative method.

The simplicity of the gradient update and its sensitivity to the condition number of the Hessian matrix are central practical features.^[ref]

A Givens rotation eliminates one component of a vector by rotating just two coordinates, preserving its Euclidean norm. For example,The length remains five: elimination moves the component into the retained coordinate rather than discarding it. Applied to a matrix, the same operation mixes two rows; QR decomposition by Givens rotations uses these local operations to build an upper triangular matrix.

For indices $i<j$, the matrix $G_{ij}(c,s)$ is the identity matrix except for its block on coordinates $i,j$:This article uses left matrix multiplication on column vectors. If $c=\cos \theta$ and $s=\sin \theta$, the block is a rotation matrix through angle $-\theta$: the minus sign belongs in the lower left. Its transpose reverses the rotation. Direct multiplication gives $G^{T}G=I$ and determinant one, so it is an orthogonal matrix; consequently $\Vert Gx{\Vert }_2^2=x^T{G}^{T}Gx=\Vert x{\Vert }_2^2$. Only two coordinates change, but this Euclidean norm identity holds for the entire vector.

The Ford-Fulkerson algorithm repeatedly augments a feasible flow along an augmenting path. When no such path remains, the source-reachable vertices in the residual network define a cut of capacity equal to the flow value, proving optimality by the max-flow min-cut theorem.

Equally spaced interpolation constructs an interpolant from data at nodes $x_j=a+jh$. It is convenient when measurements or a computational grid already have uniform spacing. The degree-at-most-$n$ polynomial interpolant through $n+1$ distinct nodes is unique, but uniqueness does not imply good approximation as $n$ increases.

For $f(x)=1/(1+25x^2)$ on $[-1,1]$, high-degree polynomial interpolation at equally spaced nodes exhibits Runge's phenomenon near the endpoints. Chebyshev nodes reduce that endpoint instability by clustering there; piecewise spline interpolation offers another alternative. The data locations and the degree must therefore be treated as approximation choices, rather than assuming that more uniformly spaced data always improve a global polynomial.

Constrained optimization minimizes or maximizes an objective function over a feasible set defined by equality or inequality constraints. In $\min f(x)$ subject to $g_i(x)\le 0$ and $h_j(x)=0$, a feasible optimum must balance improvement of the objective against the permitted directions of motion.

For example, minimizing $x^2+y^2$ subject to $x+y=1$ gives $x=y=1/2$, since $x^2+y^2=((x+y{)}^2+(x-y{)}^2)/2\ge 1/2$. The unrestricted minimizer $(0,0)$ is infeasible. Lagrange multipliers and the Karush-Kuhn-Tucker conditions generalize this balancing argument; necessity requires appropriate hypotheses, and sufficiency needs extra structure such as convexity.^[ref]

A backward differentiation formula is an implicit linear multistep method obtained by interpolating recent solution values and matching the interpolant's derivative to the ordinary differential equation at the newest time. At constant step $h$, the first-order formula is Backward Euler method; the second-order formula isThe new value appears inside $f$, so a nonlinear system may have to be solved at every step. These formulas are useful for stiff equations: their stability regions permit much larger negative-real test parameters than typical explicit methods. The BDF2 method is A-stable; this property does not extend to every higher-order formula.^[ref]