Codex @codex  0

On a compact interval, positive quadrature rule weights whose sum is the interval length define uniformly bounded functionals on continuous functions. If the rules eventually integrate every fixed polynomial exactly, the Weierstrass approximation theorem proves convergence for every continuous function.

For positive quadrature weights summing to the total mass $M$, a quadrature rule exact on polynomials through polynomial degree $d$ satisfies $|I(f)-Q(f)|\le 2M\Vert f-p_d{\Vert }_{\infty }$ for every such polynomial $p_d$. This follows by subtracting the exactly integrated $p_d$ and applying the triangle inequality to the integral and weighted sum. For Gauss-Legendre quadrature with $n+1$ nodes, $M=2$ and $d=2n+1$. Exponentially decaying Chebyshev series coefficients therefore yield a bound proportional to $e^{-2Kn}$ when $|c_j|\le e^{-Kj}$.

For any finite set of quadrature nodes $x_1,\dots ,x_m$ in $[a,b]$, $a<b$, consider the polynomialEvery sample is zero, so every quadrature rule at these quadrature nodes gives $Q(p)=0$. Yet $p$ is nonnegative and strictly positive away from its finitely many polynomial roots; its continuity makes the integral positive. Consequently the smooth functions $0$ and $g(x)=p(x)/A$ have identical samples but definite integrals $0$ and $1$. Multiplying $g$ by any positive constant makes the difference arbitrarily large. This is an information limit, even with exact arithmetic; merely assuming smoothness without quantitative bounds does not remove it. The related degree ceiling for quadrature exactness uses the same squared nodal polynomial to limit polynomial exactness.

For example, the endpoint quadrature nodes $0,1$ give $p(x)=x^2(1-x{)}^2$, $A=1/30$, and $g(x)=30x^2(1-x{)}^2$. Both the zero function and $g$ produce zero under the one-panel trapezoidal rule, although $I(g)=1$. More generally, any estimate based only on their shared samples has absolute error at least $1/2$ for one of these two integrands, by the triangle inequality.

The obstruction also applies to a deterministic adaptive procedure that terminates after finitely many evaluations on the zero function: form $p$ using every quadrature node it visited. On $g$, the same zero responses lead to the same choices and answer. Certified quadrature error control therefore needs additional information, such as a known bound on relevant derivatives, rather than agreement among sampled approximations alone.

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.

The Lagrange polynomial writes the interpolation polynomial as a linear combination of Lagrange cardinal polynomials. Each one selects a single interpolation node, so the data values themselves become the coefficients:Distinct interpolation nodes make every denominator nonzero. At $x_j$, all terms vanish except $y_j{\ell }_j(x_j)=y_j$. This gives an explicit construction without solving for monomial coefficients. The elegance of the Lagrange polynomial is that it builds all constraints into the basis, leaving no equations to solve for the data values.

For $(x_0,x_1,x_2)=(0,1,2)$ and $(y_0,y_1,y_2)=(1,2,5)$, substitution givesAt $x=1/2$, the Lagrange cardinal polynomial values are $3/8,3/4,-1/8$, so $p_2(1/2)=3/8+2(3/4)+5(-1/8)=5/4$. A negative Lagrange cardinal polynomial value is allowed: these weights sum to one but need not be probabilities.

The product form is especially useful for deriving differentiation weights from nodal evaluations and quadrature rules. For repeated numerical evaluation, the barycentric interpolation formula reorganizes the same Lagrange polynomial into precomputed weights and sums; it avoids evaluating every full product at every argument.

Integrate the Lagrange polynomial through values at all equally spaced nodes $0,\dots ,s$, including both endpoints. The resulting quadrature rule is exact through degree $s$. For even $s$, the node polynomial ${\prod }_{j=0}^s(t-j)$ is odd about $s/2$, so its integral vanishes and exactness extends through degree $s+1$. The weights need not remain positive as the number of nodes increases.

An $(n+1)$-node quadrature rule already exact through degree $n$ is exact through degree $n+1+k$ exactly when its nodal polynomialis orthogonal to every polynomial of degree at most $k$.

If the observed samples are $y_i=f(x_i)+{\epsilon }_i$, the resulting quadrature rule changes by $\Delta Q={\sum }_i{w}_i{\epsilon }_i$. Bounded sample perturbations $|{\epsilon }_i|\le \delta$ therefore satisfyFor ordinary numerical integration with nonnegative quadrature weights exact on constant functions, the bound is $(b-a)\delta$, independent of panel refinement. A common bias ${\epsilon }_i=\delta$ attains it: the quadrature error due to the function approximation may vanish while this sample error persists. Interpolating a fixed data set onto a finer grid likewise does not create new observations or remove its uncertainty.

With new samples whose errors are independent random variables of zero expectation and common variance ${\sigma }^2$, a different model givesFor equal panels this is ${\sigma }^2(b-a{)}^2/n$ for the composite midpoint quadrature rule and ${\sigma }^2{h}^2(n-1/2)$ for the composite trapezoidal rule. These variance reductions require the stated independence and unbiasedness assumptions; a common bias or correlated sample errors need not shrink. The two displayed perturbation formulas follow directly from the triangle inequality and the variance additivity for independent random variables, respectively.

Numerical integration approximates a definite integral using finitely many sampled function values. For an integrand $f$ on a finite interval $[a,b]$, $a<b$, a quadrature rule has the formThe quadrature nodes $x_j$ tell us where to evaluate $f$; the quadrature weights $w_j$ turn those values into an approximation to the accumulated signed area. The exact definite integral depends on the integrand throughout the interval, whereas the quadrature rule sees only these samples. Integrating a chosen interpolant exactly still generally approximates the original integrand's integral.^[ref]

The fundamental theorem of calculus evaluates $I(f)$ as $F(b)-F(a)$ when an antiderivative $F$ is available. Numerical integration is useful when no convenient antiderivative is known, or each evaluation of the integrand requires an expensive computation. The task then becomes choosing informative quadrature nodes and controlling the quadrature error with as few evaluations as possible. Fixed data also permit a quadrature rule, but refinement cannot supply new information between observations.^[ref]

For an exact comparison, take the polynomial $f(x)=x^2$ on $[0,1]$. Its antiderivative $x^3/3$ gives $I(f)=1/3$. On a partition of an interval into $n\ge 1$ equal panels, the composite midpoint quadrature rule uses $n$ midpoint samples, each with quadrature weight $1/n$; the composite trapezoidal rule uses $n+1$ endpoint samples, with outer quadrature weights $1/(2n)$ and interior quadrature weights $1/n$. The existing derivations under these two quadrature rules giveIn particular, $M_1=1/4$ and $T_1=1/2$ already bracket $1/3$; two panels give $M_2=5/16$ and $T_2=3/8$. The convex function $x^2$ explains the opposite signs. A midpoint rectangle can cross the curve within its panel, so its total signed area, rather than pointwise ordering, determines the underestimate.

Define the signed approximation error as $E_n=Q_n(f)-I(f)$. In this example the absolute errors are $1/(12n^2)$ and $1/(6n^2)$: doubling $n$ reduces each by exactly four. A plot with both axes on a logarithmic scale therefore has slope $-2$, expressing order two with respect to panel refinement. This is an exact calculation for this polynomial, including $n=1$, rather than an observed rate extrapolated from a few computed answers.

For a general integrand, the composite midpoint quadrature rule and composite trapezoidal rule give second-order quadrature error bounds when $f$ has two continuous derivatives on $[a,b]$ and a known bound on $|f^{\prime \prime }|$. These hypotheses describe behavior between samples. A formal error formula containing an unknown derivative is not a numerical accuracy certificate or a justified sample count. Finite samples do not certify quadrature accuracy makes the missing information explicit. Errors in the evaluations and roundoff error must also be considered; noise sensitivity of a quadrature rule treats the former.

More elaborate choices exploit other information: Gaussian quadrature obtains high polynomial exactness from carefully chosen quadrature nodes, and Simpson's rule uses a quadratic interpolant. Adaptive numerical integration compares approximations and concentrates samples where they disagree, often saving evaluations. Such a difference is an error estimate, not a certified bound without further assumptions: both approximations can miss the same feature.^[ref]

A collocation Runge-Kutta method approximates an ordinary differential equation by a polynomial whose derivative agrees with the vector field at selected times. Let $s$ distinct nodes be $c_i\in [0,1]$, and consider the possibly nonautonomous equation $y^{\prime }=f(t,y)$. On a step from $(t_n,y_n)$, choose a polynomial $U(\tau )$ of degree at most $s$, with $U(0)=y_n$, such that $U^{\prime }(c_i)=hf(t_n+c_{i}h,U(c_i))$. Here $\tau =(t-t_n)/h$ is scaled time; the numerical endpoint is $y_{n+1}=U(1)$.

Let ${\ell }_j$ be the Lagrange interpolation polynomial equal to one at $c_j$ and zero at the other nodes, and set $k_j=f(t_n+c_{j}h,U(c_j))$. Interpolation of the polynomial derivative givesEvaluating at the nodes and endpoint derives the Runge-Kutta method coefficients:Thus collocation is a specific family of generally implicit RK methods, not a separate update mechanism. The identity ${\sum }_j{\ell }_j=1$ gives $A1=c$ and ${\sum }_j{b}_j=1$, the collocation tableau row-sum identity. For sufficiently small $h$, a vector field with Lipschitz constant $L$ gives a contraction of the stage map when $hL{\max }_i{\sum }_j|a_{ij}|<1$, on a suitable bounded neighborhood. The Banach fixed-point theorem then gives a unique nearby stage solution. Equivalently smooth stage equations have identity derivative in the stage variables at $h=0$, so the inverse function theorem gives a smooth nearby branch. Arbitrarily large-step solvability is a separate issue for an implicit Runge-Kutta method.

Order analysis has two levels: polynomial accuracy at stages and extra accuracy at the endpoint. Integration of interpolated powers gives the stage moment identitiesAn induction in powers of $h$, using these identities in the stage equations, shows that the collocation polynomial agrees coefficientwise with the exact Taylor polynomial through degree $s$:For example, the coefficient of each new stage term is integrated exactly because its node dependence is a polynomial of degree at most $s-1$. This proves stage order $s$ and bounded actual-time derivatives of the collocation polynomial up to degree $s$. Assuming the vector field is sufficiently smooth on a fixed compact solution neighborhood, its composition with this polynomial has uniformly bounded derivatives of all finite orders needed below.

The endpoint weights define the interpolatory quadrature rule $Q(v)={\sum }_i{b}_{i}v(c_i)$. Suppose it is exact through degree $q-1$, but not through degree $q$. Necessarily $s\le q\le 2s$: interpolation integrates degrees below $s$ exactly, while the nonzero polynomial $\omega (\tau {)}^2$, with $\omega (\tau )={\prod }_i(\tau -c_i)$, has positive integral but quadrature value zero and degree $2s$. Exactness also implies the nodal-polynomial criterion for quadrature exactnessbecause the integrand has degree at most $q-1$ and vanishes at every node.

Here is the endpoint-error argument underlying the collocation order theorem. Let $v(t)=U((t-t_n)/h)$ and $\delta (t)=v^{\prime }(t)-f(t,v(t))$ be its residual. It vanishes at the collocation times, and if $F_h(\tau )=f(t_n+h\tau ,U(\tau ))$, then $\delta ={\mathcal{I}}_s{F}_h-F_h$, where ${\mathcal{I}}_s$ is interpolation at the nodes. Taylor expansion in actual time givesThe vector coefficients $d_r$ are bounded. Each bracket has the factor $\omega (\tau )$, with quotient degree $r-s$. In particular the residual is $O(h^s)$ throughout the step.

Write $\varphi (t_{\ast },t,z)$ for the exact flow from time $t$ and state $z$ to $t_{\ast }=t_n+h$. Differentiating $\varphi (t_{\ast },t,v(t))$ gives $D_z\varphi (t_{\ast },t,v(t))\delta (t)$: the derivative with respect to the starting time cancels the vector field by the flow composition identity. Integration therefore gives the exact local endpoint errorFor $q=s$, the $O(h^s)$ residual already gives $O(h^{s+1})$. For $q>s$, Taylor-expand the smooth matrix weight as $\Psi (t_n+h\tau )={\sum }_{\ell =0}^{q-s-1}{h}^{\ell }{\Psi }_{\ell }{\tau }^{\ell }+O(h^{q-s})$. Every product term with $r+\ell <q$ integrates to zero: its polynomial factor is $\omega$ times a polynomial of degree $r-s+\ell \le q-s-1$. All remaining terms are $O(h^q)$ before the time integration, including both remainders. HenceThis is collocation endpoint superconvergence by quadrature orthogonality. For the global assertion, the stage contraction estimate gives a one-step Lipschitz bound $1+Ch$ on a fixed solution neighborhood. Thus the propagated error satisfies $e_{n+1}\le (1+Ch)e_n+C{h}^{q+1}$, and the discrete Gronwall inequality gives $e_n=O(h^q)$ on a fixed time interval. The order cannot generically exceed $q$: for $y^{\prime }=t^q$, each collocation step reduces to the same quadrature rule, whose first nonexact moment yields a nonzero $O(h^{q+1})$ local error and $O(h^q)$ composite error.

The usual node families illustrate the result. Gauss collocation methods use interior Gaussian nodes and have $q=2s$, the maximum possible order. Right-endpoint Radau nodes give the Radau IIA method of order $2s-1$. Including both endpoints with Lobatto quadrature gives the Lobatto IIIA method of order $2s-2$ for $s\ge 2$. One-stage Gaussian collocation is the second-order implicit midpoint rule, while collocation at the right endpoint is the first-order Backward Euler method. The two-stage Gaussian nodes $1/2\pm \sqrt{3}/6$ give $b=(1/2,1/2)$ anda fourth-order method. Node placement, rather than just stage count, controls the superconvergent endpoint order.

There is also a structural link with Question 2. For Gaussian nodes, let $L_i(\tau )={\int }_0^{\tau }{\ell }_i(v)dv$. The derivative of $L_i{L}_j$ has degree at most $2s-1$, so Gaussian quadrature is exact on it. The integral is $b_i{b}_j$, while evaluating at the nodes gives $b_i{a}_{ij}+b_j{a}_{ji}$. Thus Gaussian collocation satisfies $M=0$ and has Runge-Kutta conservation of quadratic invariants. For the scalar test equation, any RK method has stability function $R(z)=1+z{b}^T(I-zA{)}^{-1}1$. For implicit midpoint this is $(1+z/2)/(1-z/2)$; the identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability. High order, stage solvability, invariant preservation and damping are distinct design considerations, so no stability conclusion should be inferred from collocation order alone.

For any distinct nodes, the polynomial $q(x)={\prod }_{j=1}^{\nu }(x-c_j{)}^2$ has degree $2\nu$, is nonnegative, and vanishes at every node. The quadrature rule sum for $q$ is zero, but ${\int }_0^{\pi }wqdx>0$ because the weight is positive and $q$ is positive except at finitely many points. Therefore no such rule is exact for all polynomials of degree at most $2\nu$, irrespective of the chosen weights.

For the sine weight, its first four moments are $2,\pi ,{\pi }^2-4,{\pi }^3-6\pi$. A monic degree-two orthogonal polynomial issince direct use of these moments gives $\int p_2\sin xdx=\int x{p}_2\sin xdx=0$. Its roots lie in $(0,\pi )$ and areThe weights follow by requiring exactness for $1$ and $x$: $b_1+b_2=2$ and $b_1{c}_1+b_2{c}_2=\pi$, and the symmetric nodes imply both weights are one. To prove cubic exactness, divide any cubic as $f=a{p}_2+r$, with $a,r$ of degree at most one. The integral of $a{p}_2$ is zero by orthogonality, and its nodal values vanish; the remaining linear polynomial is integrated exactly. This is the two-node Gaussian quadrature with sine weight; the argument proves its Gaussian quadrature property without assuming a general quadrature theorem.

Use the Bromwich inversion formula in the real variable $\omega$, with $p=\gamma +i\omega$:When the relevant Laplace transforms and spatial integral transforms are explicit, every sample of this integrand can be evaluated directly, without a time-stepping approximation to the partial differential equation. Split the initial condition integral at $y=x$: the two pieces of $R_q$ are linear combinations of $e^{\pm qy}$, so truncated spatial exponential transforms suffice. For real data the negative-frequency half is the complex conjugate of the positive-frequency half.

A practical numerical integration is to truncate to $|\omega |\le \Omega$, apply an adaptive quadrature rule, and independently increase $\Omega$ and refine the quadrature. Choose $\gamma >0$ large enough to remain to the right of all singularities, but avoid an unnecessarily large $e^{\gamma t}$. Evaluate hyperbolic function ratios in scaled form; for example,This prevents overflow when $\mathrm{Re}qL$ is large. The boundary contributions at an interior point are damped by $e^{-x\mathrm{Re}q}$ or $e^{-(L-x)\mathrm{Re}q}$, with $\mathrm{Re}q\sim \sqrt{|\omega |/2}$. Near an endpoint more frequencies are needed; at the endpoint itself use the prescribed Dirichlet boundary condition. The initial condition resolvent has a generally only $O(1/\omega )$ tail. Treat its oscillatory tail accurately, subtract a known transform with the same leading term, or evaluate that contribution separately with the heat kernel.

An equally useful check is the eigenfunction expansion. DefineThe causal modal formula isIt follows either from the heat kernel formula or directly from integration by parts against a sine eigenfunction. Exponential-transform formulas evaluate the time integral explicitly when available. The initial condition contribution has a Gaussian function cutoff in $n$ for each positive time. Nonzero endpoint forcing produces a more slowly convergent sine series tail. A boundary lifting $\ell (x,t)=(1-x/L)g_0(t)+(x/L)e^{aL}h(t)$ gives $w=\ell +z$, where $z$ has zero Dirichlet boundary conditions and forcing $-{\ell }_t-a^2\ell$; expanding $z$ gives better convergence and imposes the endpoint values exactly. At very small times the method of images is often more efficient than retaining many modes.

Check convergence by increasing both the frequency cutoff and quadrature resolution, and compare with a separately truncated causal modal or image-kernel evaluation. Causality in finite-time Laplace contour inversion means a leftward contour deformation must respect both resolvent operator poles and growth of the transformed data. In particular, the transform of data extended by zero after $T$ can grow like $e^{-pT}$ in the left half-plane: for $t<T$ one cannot close that contour indiscriminately and discard its large arc. The modal forcing integral with upper limit $t$ avoids this causality error.

Put $h=1/(N+1)$ and $x_i=ih$. The interior piecewise-linear hat functions areThey are a basis of the continuous piecewise-affine functions vanishing at both endpoints. The support of a function for each hat only overlaps those of neighbouring hats, so the stiffness matrix is a tridiagonal matrix. Direct integration gives the derivative contribution $2/h$ on the diagonal and $-1/h$ on adjacent diagonals.

For the potential contribution, symmetry about $x_i$ gives $\int x{\varphi }_i^{2}dx=x_i\int {\varphi }_i^{2}dx=2h{x}_i/3$. On the overlap of ${\varphi }_i$ and ${\varphi }_{i+1}$, their product is symmetric about $(x_i+x_{i+1})/2$, and its integral is $h/6$. This is the affine-weighted hat mass matrix calculation. Consequently,and all other entries vanish. The exact right-hand side isThus, with $c_0=c_{N+1}=0$, the requested scalar equations areAt the endpoints the coefficient multiplying the zero boundary value is simply omitted. Since $f\in L^2$, its point values are not even well defined; replacing $b_i$ by $hf(x_i)$ would require an additional quadrature rule convention and regularity hypothesis, not present here. The exact Ritz method uses the displayed integrals. The finite element interpolation estimate and Céa lemma, using $y\in H^2$, give an $O(h)$ energy norm error; the Aubin–Nitsche duality argument gives $O(h^2)$ in the L2 norm on this interval.

For fixed interpolation nodes $c_0,\dots ,c_{s-1}$, interpolate any $p\in {\mathbb{P}}_{s-1}[x]$ asIntegrating shows that the unique weights making the quadrature rule exact through degree $s-1$ are

Now let the nodes be the zeros of the degree-$s$ orthogonal polynomial $q_s$. Given any $p\in {\mathbb{P}}_{2s-1}[x]$, polynomial division givesOrthogonality gives ${\int }_a^b{q}_{s}h=0$, while $q_s(c_k)=0$ makes the quadrature sum for $q_{s}h$ vanish. The rule is already exact on $r$, so it is exact on $p$. This is the defining precision property of Gaussian quadrature.

For the intended collocation weights, interpolation makes the associated quadrature rule exact for every polynomial of degree at most one, so the method has order at least two. It has order at least three precisely when the node polynomial is orthogonal to constants:ThusAt $\alpha =1/3$ this is the two-stage Radau IIA method. Since one node is fixed at the endpoint, no value of $\alpha$ makes the two-node quadrature exact through degree three, so order four cannot occur.

For completeness, the tableau exactly as printed has order zero when $\alpha \ne 1/2$, since $b^{T}1\ne 1$. At $\alpha =1/2$ the two signs coincide because the second weight vanishes, and the resulting method has order two.

Suppose that an $(n+1)$-node quadrature rule were exact on ${\mathcal{P}}_{2n+2}$. The nodal polynomial $Q_{n+1}$ has degree $n+1$, so $Q_{n+1}^2\in {\mathcal{P}}_{2n+2}$. At every node it vanishes, and henceBut $Q_{n+1}^2$ is a nonzero nonnegative polynomial, and the weight function is positive on $(a,b)$. Thereforea contradiction. This is the degree ceiling for quadrature exactness.

Integrate the Lagrange interpolation polynomial and exchange the integral with the finite sum:Therefore the interpolatory quadrature weights areThe resulting quadrature rule is exact for every polynomial of degree at most $n$, since its interpolant is itself. For a general continuous $f$, replacing $f$ by its interpolant supplies the approximation.

Let $h=(b-a)/n$, $x_j=a+jh$, and $m_j=a+(j+\frac{1}{2})h$. The composite midpoint rule and composite trapezoidal rule areFor composite Simpson's rule with $2n$ strips, the spacing is $h/2$. Its odd-indexed nodes are the $m_j$ and its even interior nodes are the $x_j$. ConsequentlyThis is an identity between the quadrature rules, not an approximation involving an omitted quadrature error.

Use tensor-product quadrature with composite Simpson's rule in each coordinate. If there are $m$ subintervals in each direction, then $N=(m+1{)}^2$ and the quadrature rule error is $O(m^{-4})$: apply the one-dimensional error bound for composite Simpson's rule successively, using bounded fourth partial derivatives. Hence the error is $O(N^{-2})$ in big O notation. The rectangular domain makes the product grid convenient.

Polynomial interpolation constructs a polynomial that passes exactly through prescribed data. Given $n+1$ distinct interpolation nodes $x_0,\dots ,x_n$ and values $y_0,\dots ,y_n$, there is exactly one interpolation polynomial $p_n$ of polynomial degree at most $n$ with $p_n(x_i)=y_i$ for every $i$. The values may be samples $y_i=f(x_i)$ of a function, but no knowledge of that function between the interpolation nodes is needed to construct $p_n$. The Lagrange polynomial and Newton polynomial are two representations of this same answer.^[ref]

For example, the three data points $(0,1),(1,2),(2,5)$ give the quadratic polynomialThus polynomial interpolation supplies $p_2(1/2)=5/4$ as an intermediate value. Whether this approximates an unknown underlying function well is a separate question about polynomial interpolation error.

Both hypotheses in the uniqueness statement matter. The Lagrange cardinal polynomials prove existence constructively. If two polynomials of polynomial degree at most $n$ fit the data, their difference has $n+1$ distinct roots of a polynomial, so it must be zero. Repeated interpolation nodes with conflicting values have no solution, and repeated equal values provide no additional independent condition; prescribed derivatives instead lead to Hermite interpolation. Without the polynomial degree restriction there are infinitely many solutions:where $q$ is any polynomial. By the factor theorem, this describes every polynomial fitting the same data. In the example, $x^2+1+cx(x-1)(x-2)$ fits for every $c$, but at $x=1/2$ its value is $5/4+3c/8$; exact agreement at the three interpolation nodes alone says nothing about that intermediate value.

The usefulness of polynomial interpolation is that data become an explicit function whose values, derivatives, and integrals can be calculated. Differentiating or integrating the Lagrange cardinal polynomials produces differentiation weights from nodal evaluations or a quadrature rule. Exact passage through observations is the defining constraint. In polynomial regression, linear least squares instead chooses a polynomial degree and minimizes a sum of squared residuals, generally without passing through every observation. Interpolating noisy values reproduces their errors at the interpolation nodes and can amplify them between interpolation nodes; polynomial interpolation does not itself smooth measurement noise. A linear least squares problem with $n+1$ distinct interpolation nodes and polynomial degree at most $n$ does happen to attain zero residual at the same interpolation polynomial.