A Newton-Cotes quadrature integrates a Lagrange interpolation polynomial at equally spaced quadrature nodes. For the Newton-Cotes closed quadrature, the endpoints are included: write $x_j=a+jh$, $h=(b-a)/n$, and useThese quadrature weights depend on the quadrature nodes, not on the sampled function. Open Newton-Cotes quadrature uses equally spaced interior quadrature nodes instead.

For Simpson's rule, use $n=2$ and set $x=a+ht$, so $0\le t\le 2$. The Lagrange cardinal polynomials becomeIntegrating these quadratic polynomials gives $1/3$, $4/3$, $1/3$, respectively. The change of variable contributes $h$, soIts degree of precision of a quadrature rule is three: it integrates every quadratic polynomial by construction, and also the cubic nodal polynomial $(x-a)(x-a-h)(x-b)$, whose integral vanishes by odd symmetry about $a+h$. Any cubic polynomial differs from its quadratic polynomial interpolant by a constant multiple of this nodal polynomial.

The product rule applied at the interpolation node $x_k$ givesThus Lagrange interpolation gives, away from the interpolation nodes,Apply the same Lagrange interpolation formula to the constant polynomial one. Since its interpolant is itself,Dividing the preceding two identities proves the barycentric interpolation formulaThe denominator is exactly $1/w(x)$, so is nonzero away from the interpolation nodes. At $x=x_j$, multiply numerator and denominator by $x-x_j$ and take the limit: their values tend to $w_j{f}_j$ and $w_j$, respectively. The removable singularity therefore has value $f_j$; in numerical evaluation return the stored value directly at an interpolation node.

For each interpolation node define the Lagrange cardinal polynomialAll denominators are nonzero because the interpolation nodes are distinct. Each polynomial has polynomial degree at most $n$, and ${\ell }_j(x_i)={\delta }_{ij}$, where ${\delta }_{ij}$ is the Kronecker delta. Consequently the linear combinationbelongs to the vector space of univariate polynomials ${\Pi }_n$ and has the required values. This proves existence of the Lagrange interpolation polynomial.

For uniqueness, suppose $p_n$ and $r_n$ both work. Their difference is a polynomial of polynomial degree at most $n$ with $n+1$ distinct roots of a polynomial. The root bound for a polynomial forces $p_n-r_n$ to be the zero polynomial.

The seven supplied values follow $f_i=i^2+1$. The quadratic polynomial $x^2+1$ belongs to ${\Pi }_6$ and takes all seven values, so uniqueness givesThe subscript bounds the polynomial degree; it does not require polynomial degree exactly six.

Since $Q^{-1}=Q^T$, the QR factorization givesFirst compute $y=Q^{T}b$ by a matrix-vector product, costing $n^2$ scalar multiplications and $n(n-1)$ additions for dense $Q$. Then solve the upper triangular matrix system $Rx=y$ by back substitution:Every denominator is nonzero because $R=Q^{T}A$ is invertible. The back substitution costs $n(n-1)$ multiplications and subtractions together, plus $n$ divisions, for $n^2$ arithmetic operations. Including $Q^{T}b$, the total is $3n^2-n$ arithmetic operations under this scalar counting convention, hence $O(n^2)$, with $O$ denoting Big O notation. The already supplied QR factorization is reused; no matrix inverse needs to be computed.

A real square matrix $Q$ is an orthogonal matrix when $Q^{T}Q=I$, equivalently its columns form an orthonormal basis. Here $Q^T$ is the transpose and $I$ the identity matrix. For every vector $x$,Taking nonnegative square roots provesThus an orthogonal matrix preserves the Euclidean norm.

Assume $f\in C^{n+1}[a,b]$ and set $\omega (t)={\prod }_{i=0}^n(t-x_i)$, a monic polynomial of polynomial degree $n+1$. Fix $x\in [a,b]$ distinct from all the interpolation nodes, and chooseThe auxiliary functionvanishes at the $n+2$ distinct points $x_0,\dots ,x_n,x$. Order those points and repeatedly apply Rolle's theorem. After $n+1$ applications, there is some $\xi$ strictly between their smallest and largest values, hence in $(a,b)$, such that $F^{(n+1)}(\xi )=0$.

Because $p_n$ has polynomial degree at most $n$ and $\omega$ is monic of polynomial degree $n+1$,Substituting the resulting value of $C$ proves the polynomial interpolation error formulaAt an interpolation node, both sides are zero, so any $\xi \in (a,b)$ works. Such an interior point exists because $n\ge 1$ and the interpolation nodes are distinct. The selected $\xi$ generally depends on $x$; the theorem makes an existence assertion, not a regularity assertion about that dependence.

The polynomial $p_n$ interpolates $f$ at the interpolation nodes whenFor distinct interpolation nodes, define the Lagrange cardinal polynomialsThe denominators are nonzero, each polynomial degree is at most $n$, and ${\ell }_i(x_j)={\delta }_{ij}$, with ${\delta }_{ij}$ the Kronecker delta. Thereforeis a Lagrange interpolation polynomial belonging to the vector space of univariate polynomials ${\Pi }_n$, and it has all the required sampled values.

If $r_n\in {\Pi }_n$ also works, $p_n-r_n$ has $n+1$ distinct roots of a polynomial while its polynomial degree is at most $n$. The root bound for a polynomial forces it to be the zero polynomial, proving uniqueness. The phrase “degree $n$” denotes the allowed maximum; the actual polynomial degree may be smaller. No continuity assumption is needed merely to construct the interpolation polynomial from these finitely many values.

Write $E(h)$ for the supplied approximation. The asymptotic expansion iswhere $O$ is Big O notation. At the refined parameter $h/3$,Multiply this identity by nine and subtract the first one. The unknown leading truncation error coefficient cancels:Thus the Richardson extrapolation estimate isThe two supplied values giveThe correction to the finer approximation is $(3-2)/8=1/8$. Neglecting the higher-order remainder would also give $c=9/8$.

The conclusion is asymptotic: the leading error order improves from two to at least three when both parameter values are in the range of the stated asymptotic expansion. With no bound on the remainder constant, the two finite observations alone do not guarantee that $25/8$ is closer to $Y$. Indeed, for any proposed value $L$, the polynomial modelsatisfies $E(1)=2$, $E(1/3)=3$ and the required asymptotic expansion with true limit $L$. This is why finite observations do not certify Richardson extrapolation accuracy: asymptotic cancellation provides an improved-order estimator, while a numerical error certificate needs quantitative remainder information.