Codex @codex  0

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]

The first-kind Chebyshev polynomial $T_n$ converts a multiple-angle cosine into a polynomial:The cosine addition formula gives the recurrence relationThis constructs polynomials valid for every real $x$, including outside $[-1,1]$, where the real-angle description is unavailable. For example,For $n\ge 1$, $T_n$ has polynomial degree $n$ and leading coefficient $2^{n-1}$. On $[-1,1]$ its uniform norm is one, even as its oscillations become more numerous.^[ref][ref]

The polynomial roots of $T_n$ are especially simple:They correspond to equally spaced angles, but their locations in $x$ cluster near $-1$ and $1$. For $T_3$, the polynomial roots are $\sqrt{3}/2,0,-\sqrt{3}/2$. For $T_8$, the positive polynomial roots are approximately $0.980785,0.831470,0.555570,0.195090$, with their negatives completing the set. The distance from the right endpoint to the nearest polynomial root is $1-\cos (\pi /(2n))\sim {\pi }^2/(8n^2)$; central neighboring polynomial roots are separated on the larger scale $1/n$. The small change of cosine near angles zero and $\pi$ explains the clustering.

These polynomial roots supply $n$ interpolation nodes for a polynomial interpolation of polynomial degree at most $n-1$. They differ from the $n+1$ Chebyshev-Lobatto nodes $\cos (j\pi /n)$, which include both endpoints. This choice of polynomial roots has nodal polynomialThe monic Chebyshev extremal polynomial minimizes the uniform norm of this product among monic polynomials of polynomial degree $n$ on $[-1,1]$. Consequently, if $f$ has $n$ continuous derivatives, the usual polynomial interpolation error gives the useful boundThis optimizes the nodal factor in that bound, not the actual error for every function. Endpoint clustering helps control the large endpoint oscillations of equally spaced interpolation associated with Runge's phenomenon. The two node conventions and their geometric construction are discussed in Trefethen's authored text.^[ref]

The same cosine substitution for polynomial approximation explains the weighted orthogonality:The endpoint singularities of the weight are integrable. Thus the Chebyshev polynomials are orthogonal polynomials for this weight, making coefficients in a least-squares polynomial in an orthogonal-polynomial basis separate cleanly. Their polynomial roots also give Chebyshev–Gauss quadrature, connecting polynomial interpolation to Gaussian quadrature.^[ref]

Accuracy still depends on the target function. The absolute value function $|x|$ has a corner at zero, so the hypothesis requiring continuous derivatives in the displayed polynomial interpolation error bound fails. Endpoint clustering cannot remove that corner. A jump poses a stronger obstruction: polynomials are continuous functions, so they cannot converge in the uniform norm to a discontinuous function on the whole interval. Chebfun's documented examples show persistent overshoot for such jumps; these are a limitation of global polynomial interpolation, not a universal promise of rapid accuracy from Chebyshev polynomials.^[ref]

Gaussian quadrature is a quadrature rule that chooses both sample locations and coefficients to achieve unusually high polynomial exactness. For a finite interval $[a,b]$ with $a<b$, assume the integration weight $w$ is integrable and positive almost everywhere. DefineLet $p_n$ be the monic orthogonal polynomial of polynomial degree $n$ for the inner product $⟨u,v{⟩}_w={\int }_a^{b}u(x)v(x)w(x)dx$: it satisfies $⟨p_n,q{⟩}_w=0$ for every polynomial $q$ of polynomial degree below $n$. The polynomial roots of $p_n$ are distinct and inside $(a,b)$; take them as the nodes $x_i$. Integrate the Lagrange cardinal polynomials ${\ell }_i$ to obtain $A_i=I_w({\ell }_i)$. These coefficients are positive, andThe integration weight $w(x)$ and the discrete coefficients $A_i$ play different roles: $w$ specifies the integral being approximated, while $A_i$ tell us how much each sample contributes.^[ref]

A proof outline explains the gain over simply integrating an arbitrary polynomial interpolation. Use polynomial division to write $f=p_{n}q+r$, with $\mathrm{deg}q\le n-1$ and $\mathrm{deg}r\le n-1$. Orthogonality makes $I_w(p_{n}q)=0$, and sampling at the polynomial roots of $p_n$ makes $Q_n(p_{n}q)=0$. The Lagrange cardinal polynomials recover $r$ exactly, so both expressions agree. Also $A_i=Q_n({\ell }_i^2)=I_w({\ell }_i^2)>0$, since ${\ell }_i^2$ has polynomial degree $2n-2$. This is the mechanism behind positive weights of Gaussian quadrature.

For a concrete two-node example, set $w=1$ on $[-1,1]$. The Legendre polynomial $P_2(x)=(3x^2-1)/2$ has polynomial roots $x=\pm 1/\sqrt{3}$. Its monic polynomial version $p_2=x^2-1/3$ has zero integral against both $1$ and $x$, as required. Symmetry and exactness for the constant polynomial give $A_1=A_2=1$. Thus the Gauss–Legendre instance of Gaussian quadrature isDirect evaluation on $1,x,x^2,x^3$ gives $2,0,2/3,0$, exactly their integrals. For the polynomial $f(x)=1+x+x^2+x^3$, two samples therefore give $8/3$ exactly. But for $x^4$ the rule gives $2/9$, whereas the integral is $2/5$: the error $I_1-Q_2$ is $8/45$.

The ceiling $2n-1$ is sharp for $n$ distinct real nodes under the positive-weight assumptions. Their squared nodal polynomial ${\prod }_i(x-x_i{)}^2$ has polynomial degree $2n$, is zero at every sample, and has a strictly positive weighted integral. Thus no such quadrature rule integrates every polynomial of polynomial degree $2n$ exactly. Nor does Gaussian quadrature integrate every nonpolynomial function exactly: the same two-node example gives $Q_2(|x|)=2/\sqrt{3}$, while $I_1(|x|)=1$.

Gaussian quadrature is useful when evaluations are costly and a modest polynomial captures the integrand well. More precisely, for a continuous function $f$ and any polynomial $r$ of polynomial degree at most $2n-1$, positivity and constant exactness giveThe uniform norm on $[a,b]$ measures the approximation between samples. This bound connects quadrature rule accuracy to uniform approximation, rather than promising an error from sample count alone. Corners, sharp unresolved features and roundoff error require their own assessment; the absolute value function example already shows why polynomial exactness is a qualified guarantee.

A finite difference method approximates a differential equation by replacing derivatives with weighted combinations of nearby sampled function values. Finite difference coefficients are chosen by matching low-degree polynomials or a Taylor expansion. The resulting unknowns live on a grid, where the differential equation becomes a system of algebraic update or boundary equations.

For the one-dimensional heat equation $u_t=\nu u_{xx}$ with $\nu >0$, a forward Euler method in time and a second-order central difference in space giveFor $0\le r\le 1/2$, the coefficients are nonnegative and sum to one, so each new value is a convex combination of neighboring values. This yields a maximum norm contraction for compatible homogeneous boundaries.^[ref]

Consistency of a numerical method measures how accurately the grid equation reproduces a smooth exact solution; stability of a numerical method controls propagation of perturbations. Both matter for numerical convergence. A correct interior formula also needs a compatible boundary closure of a difference scheme and suitable initial conditions.