43.09.03 · numerical-analysis / 09-numerical-quadrature

Gauss quadrature via orthogonal polynomials

shipped3 tiersLean: none

Anchor (Master): Davis-Rabinowitz 1984 *Methods of Numerical Integration* 2e (Academic Press) Ch. 2-3 (Gauss rules, the orthogonal-polynomial connection, the Radau and Lobatto variants, and the Stieltjes error theory); Gautschi 2004 *Orthogonal Polynomials: Computation and Approximation* (Oxford) Ch. 1-3 (the Jacobi matrix, the Golub-Welsch eigenvalue algorithm, and the classical Gauss families); Gautschi 1981 'A survey of Gauss-Christoffel quadrature formulae' in *E. B. Christoffel* (Birkhäuser) (the comprehensive survey of Gauss-Christoffel rules, the Radau/Lobatto modifications, and the weight-positivity and convergence theory)

Intuition Beginner

When you estimate the area under a curve from a handful of sample heights, you face two separate choices: where to put the sample points, and how much weight to give each one. The classical rules from the previous unit fix the points in advance — equally spaced — and only get to choose the weights. That throws away half of your freedom. The Gauss idea is to choose both at once: pick the cleverest possible sample locations and the cleverest possible weights together, so the rule squeezes out as much accuracy as a fixed number of samples can possibly buy.

How much does that extra freedom win you? With $n$ samples you have $n$ locations and $n$ weights to set, so $2 n$ knobs in total. Each knob can be tuned to make the rule exact on one more polynomial. A rule with $n$ equally spaced points integrates polynomials up to degree about $n - 1$ exactly; a Gauss rule with the same $n$ points reaches degree $2 n - 1$ . You roughly double the accuracy for the same number of function evaluations, which is why Gauss rules are the workhorse behind most automatic integration software.

The surprise is where the best sample points turn out to sit. They are not spread evenly, and they are not chosen by trial and error. They are the roots of a special polynomial — the one built to be perpendicular, in a weighted-average sense, to all lower-degree polynomials. Those orthogonal polynomials came up when fitting curves by least squares, and the same family reappears here as the secret map of where to sample. One construction, two uses.

A second piece of good news is that the weights always come out positive. Positive weights mean the rule averages your sample heights rather than amplifying their errors, so unlike the high-point evenly-spaced rules, Gauss rules stay stable no matter how many points you use. More points always means more accuracy, never a blow-up.

Visual Beginner

Picture the same interval sampled two ways. The evenly spaced rule drops its points at fixed, regular positions. The Gauss rule pulls its points inward, away from the endpoints, clustering them where they do the most good — exactly at the roots of the special perpendicular polynomial. With the same count of points, the Gauss arrangement integrates far higher-degree curves with no error.

The table contrasts the two families by what they fix and what they win.

rule family	points	weights	highest degree exact ( $n$ points)
evenly spaced (Newton-Cotes)	fixed, equal spacing	chosen	about $n - 1$
Gauss	chosen (orthogonal-polynomial roots)	chosen	$2 n - 1$

Reading across: the Gauss rule spends its extra freedom — the choice of where to sample — and is repaid with roughly twice the exact degree. The points land at the orthogonal-polynomial roots, which always sit strictly inside the interval.

Worked example Beginner

Let us build the two-point Gauss rule on the interval from $- 1$ to $1$ with the plain weight (every part of the interval counts equally), and check it integrates a cubic perfectly.

Step 1. Find the sample points. For two points the special perpendicular polynomial is $x^{2} - 1 \div 3$ . Its roots are where $x^{2} = 1 \div 3$ , namely $x = - 1 \div 3$ and $x = + 1 \div 3$ , about $- 0.577$ and $+ 0.577$ . These are the two Gauss points — pulled inward from the endpoints.

Step 2. Find the weights. By symmetry the two weights are equal, and they must add up to the length of the interval, which is $2$ . So each weight is $1$ . The rule is: add the function value at $- 0.577$ and the function value at $+ 0.577$ .

Step 3. Test it on $f (x) = x^{3}$ . The true area under $x^{3}$ from $- 1$ to $1$ is $0$ , because the curve is odd and the two halves cancel. The rule gives $(- 0.577)^{3} + (0.577)^{3} = 0$ , since the two cubes cancel. Match.

Step 4. Test it on $f (x) = x^{2}$ . The true area is the value of $x^{3} \div 3$ at the endpoints: $1 \div 3 - (- 1 \div 3) = 2 \div 3$ . The rule gives $(- 0.577)^{2} + (0.577)^{2} = 1 \div 3 + 1 \div 3 = 2 \div 3$ . Match again.

Step 5. Confirm the badge. Two points buy exactness up to degree $2 \times 2 - 1 = 3$ . A check on $x^{3}$ and $x^{2}$ both passed, and by the same cancellation the rule is exact for every cubic.

What this tells us: with only two well-placed samples and two equal weights, the Gauss rule integrates every cubic exactly — the same feat Simpson's rule needs three points for. Choosing the points, not just the weights, is what pays for the extra degree.

Check your understanding Beginner

Formal definition Intermediate+

Throughout, $w$ is a fixed weight function on an interval $[a, b]$ (possibly infinite) with $w (x) > 0$ almost everywhere and all moments $μ_{k} = \int_{a}^{b} x^{k} w (x) d x$ finite, exactly the setting of 43.08.05. The weighted inner product is $(f, g)_{w} = \int_{a}^{b} w (x) f (x) g (x) d x$ , $Π_{k}$ is the space of real polynomials of degree at most $k$ , and $ϕ_{0}, ϕ_{1}, ϕ_{2}, \dots$ with $de g ϕ_{k} = k$ are the orthogonal polynomials for $w$ of 43.08.05, with leading coefficient $k_{n}$ in $ϕ_{n}$ . The error functional of a rule $Q$ is $E (f) = \int_{a}^{b} w f d x - Q (f)$ , as in 43.09.01. Notation used below — $\int$ , $\prod$ , $\sum$ , $\partial$ , the orthogonal polynomial $ϕ_{n}$ , the leading coefficient $k_{n}$ , and the Jacobi matrix $J_{n}$ — is recorded in _meta/NOTATION.md.

Definition (general quadrature rule and the optimal-exactness problem). An $n$ -point quadrature rule for the weight $w$ is $Q (f) = \sum_{i = 1}^{n} w_{i} f (x_{i})$ with $n$ nodes $x_{1} < \dots < x_{n}$ in $[a, b]$ and $n$ weights $w_{i}$ , approximating $\int_{a}^{b} w f d x$ . The rule has $2 n$ free parameters. The optimal-exactness problem is to choose the nodes and weights so that $E (p) = 0$ for every $p \in Π_{d}$ with $d$ as large as possible. Since $2 n$ parameters can match at most $2 n$ conditions, the ceiling is $d = 2 n - 1$ .

Definition (interpolatory weights). Given fixed nodes $x_{1}, \dots, x_{n}$ , the interpolatory weights are $w_{i} = \int_{a}^{b} w (x) L_{i} (x) d x$ , where $L_{i}$ are the Lagrange cardinal polynomials on the nodes from 43.08.01; this is the weight choice of 43.09.01 that makes the rule exact on $Π_{n - 1}$ for any node set.

Definition (Gauss quadrature rule). The Gauss (Gauss-Christoffel) quadrature rule of order $n$ for $w$ takes its nodes to be the $n$ roots of the orthogonal polynomial $ϕ_{n}$ — real, simple, and inside $(a, b)$ by 43.08.05 — and its weights to be the interpolatory weights $w_{i} = \int_{a}^{b} w L_{i} d x$ on those nodes. The weights $w_{i}$ are also called Christoffel numbers.

Definition (Gauss families). Specialising the weight gives the classical families: $w \equiv 1$ on $[- 1, 1]$ gives Gauss-Legendre; $w (x) = (1 - x^{2})^{- 1/2}$ on $[- 1, 1]$ gives Gauss-Chebyshev; $w (x) = e^{- x^{2}}$ on $(- \infty, \infty)$ gives Gauss-Hermite; $w (x) = e^{- x}$ on $[0, \infty)$ gives Gauss-Laguerre. The Gauss-Radau rule fixes one endpoint as a node and is exact on $Π_{2 n - 2}$ ; the Gauss-Lobatto rule fixes both endpoints and is exact on $Π_{2 n - 3}$ .

Counterexamples to common slips Intermediate+

"More nodes can always raise the exactness past $2 n - 1$ ." No node-and-weight choice with $n$ nodes beats $2 n - 1$ : the polynomial $\prod_{i} (x - x_{i})^{2} \in Π_{2 n}$ is strictly positive away from the nodes, so its weighted integral is positive while every rule evaluates it to $0$ . The ceiling is structural, not a matter of cleverness.
"The Gauss nodes are equally spaced, just rescaled." They are the roots of $ϕ_{n}$ , which are not equispaced — for Gauss-Legendre they cluster toward the endpoints like $cos$ of equispaced angles. Equispacing the nodes throws away the orthogonality that buys the extra $n$ degrees.
"Positivity of the weights needs the weight $w$ to be constant." The argument $w_{i} = \int w L_{i}^{2} > 0$ uses only $w > 0$ and exactness on $Π_{2 n - 2}$ , so it holds for every admissible weight, including the singular Chebyshev weight and the infinite-interval Hermite and Laguerre weights.
"The error term uses $f^{(n + 1)}$ , like interpolatory quadrature." The Gauss error involves $f^{(2 n)}$ , the doubled order, because the rule is exact through degree $2 n - 1$ . Using the interpolation-error derivative $f^{(n + 1)}$ would understate the order by $n - 1$ .

Key theorem with proof Intermediate+

The signature result is the maximal-exactness theorem: among all $n$ -point rules, the Gauss rule alone reaches degree of exactness $2 n - 1$ , and it does so precisely by placing the nodes at the roots of $ϕ_{n}$ . The proof rests on one algebraic move — dividing a high-degree polynomial by $ϕ_{n}$ — and one analytic fact — that $ϕ_{n}$ is orthogonal to everything of lower degree. This follows the development in Süli-Mayers ^{[Süli-Mayers 2003 §10.2]}.

Theorem (Gauss quadrature: maximal degree of exactness). Let the nodes $x_{1}, \dots, x_{n}$ be the roots of $ϕ_{n}$ and the weights be the interpolatory weights $w_{i} = \int_{a}^{b} w L_{i} d x$ . Then the rule $Q (f) = \sum_{i = 1}^{n} w_{i} f (x_{i})$ satisfies $E (p) = 0$ for every $p \in Π_{2 n - 1}$ . Conversely, no $n$ -point rule is exact on $Π_{2 n}$ , and a rule is exact on $Π_{2 n - 1}$ only if its nodes are the roots of $ϕ_{n}$ .

Proof. Take any $p \in Π_{2 n - 1}$ . Divide $p$ by the degree- $n$ polynomial $ϕ_{n}$ : there are unique $q, r \in Π_{n - 1}$ with $$ p = q,\phi_n + r, \qquad \deg q \le n-1, \quad \deg r \le n-1, $$ since $de g p \leq 2 n - 1$ and $ϕ_{n}$ is monic-up-to-scale of degree $n$ . Integrate against $w$ : $$ \int_a^b w,p,dx = \int_a^b w,q,\phi_n,dx + \int_a^b w,r,dx = (\phi_n, q)w + \int_a^b w,r,dx . $$ Because $q \in \Pi{n-1} $an d$ \phi_n $i sor t h o g o na l t o$ \Pi_{n-1} $(t h e d e f inin g p r o p er t y f r o m [43.08.05]),$ (\phi_n, q)w = 0 $, so$ \int_a^b w,p,dx = \int_a^b w,r,dx $. N o w a ppl y t h er u l e . S in cee a c hn o d e$ x_i $i s a r oo t o f$ \phi_n $,$ \phi_n(x_i) = 0 $, h e n ce$ p(x_i) = q(x_i)\phi_n(x_i) + r(x_i) = r(x_i)$, and $$ Q(p) = \sum_i w_i p(x_i) = \sum_i w_i r(x_i) = Q(r). $$ The rule is interpolatory on $n$ nodes, so it is exact on $\Pi{n-1} $, an d$ r \in \Pi_{n-1} $g i v es$ Q(r) = \int_a^b w,r,dx $. C hainin g t h e tw oe q u a l i t i es,$ Q(p) = \int_a^b w,r,dx = \int_a^b w,p,dx $, so$ E(p) = 0 $. A s$ p \in \Pi_{2n-1} $w a s a r bi t r a r y, t h er u l e i se x a c t o n$ \Pi_{2n-1}$.

For the ceiling, the polynomial $g (x) = \prod_{i = 1}^{n} (x - x_{i})^{2} \in Π_{2 n}$ satisfies $g (x_{i}) = 0$ for every node, so $Q (g) = 0$ ; but $g > 0$ except at the nodes, so $\int_{a}^{b} w g d x > 0$ by positivity of $w$ . Thus $E (g) \neq = 0$ and no $n$ -point rule is exact on $Π_{2 n}$ , whatever its nodes. Finally, suppose an $n$ -point rule with nodes $y_{1}, \dots, y_{n}$ is exact on $Π_{2 n - 1}$ ; let $ψ (x) = \prod_{i} (x - y_{i}) \in Π_{n}$ . For any $s \in Π_{n - 1}$ the product $s ψ \in Π_{2 n - 1}$ vanishes at every node, so $0 = Q (s ψ) = \int_{a}^{b} w s ψ d x = (ψ, s)_{w}$ . Hence $ψ$ is orthogonal to all of $Π_{n - 1}$ ; a monic degree- $n$ polynomial with that property is $ϕ_{n}$ up to scale (uniqueness of the monic orthogonal polynomial, 43.08.05), so the nodes $y_{i}$ are the roots of $ϕ_{n}$ . $□$

Corollary (positivity of the Gauss weights). The Gauss weights are strictly positive, $w_{i} > 0$ for every $i$ , and they sum to $\int_{a}^{b} w d x$ . Fix $i$ and let $L_{i}$ be the cardinal polynomial; then $L_{i}^{2} \in Π_{2 n - 2} \subset Π_{2 n - 1}$ , so the rule integrates it exactly: $$ 0 < \int_a^b w,L_i^2,dx = Q(L_i^2) = \sum_j w_j L_i(x_j)^2 = w_i , $$ the last step because $L_{i} (x_{j}) = δ_{ij}$ . Exactness on the constant $1 \in Π_{0}$ gives $\sum_{i} w_{i} = \int_{a}^{b} w d x$ . Positive weights summing to the total mass are what make Gauss rules stable: $\sum_{i} ∣ w_{i} ∣ = \sum_{i} w_{i} = \int_{a}^{b} w$ , the minimal operator norm, in sharp contrast to the unbounded $\sum ∣ w_{i} ∣$ of high-order Newton-Cotes in 43.09.01.

Bridge. The division $p = q ϕ_{n} + r$ is the foundational reason Gauss quadrature works: it splits any polynomial of degree up to $2 n - 1$ into a part that the orthogonality of $ϕ_{n}$ annihilates under the integral and a part of degree $< n$ that the interpolatory rule handles exactly. This is exactly the payoff of the orthogonal-polynomial root structure of 43.08.05: the roots are not merely real and simple, they are the unique node set that makes the integral of $q ϕ_{n}$ vanish, which is what doubles the exactness from $n - 1$ to $2 n - 1$ . The result builds toward the error term and the Golub-Welsch algorithm of the Advanced section, and it appears again in the spectral and Gaussian-process methods of 37.08.01 and the Clenshaw-Curtis comparison of 43.09.02. The central insight is that exactness on $Π_{2 n - 1}$ is equivalent to the node polynomial being orthogonal to $Π_{n - 1}$ , which generalises the interpolatory exactness of 43.09.01 from degree $n - 1$ to degree $2 n - 1$ and is dual to the least-squares projection of 43.08.05, where the same $ϕ_{n}$ is the residual direction. Putting these together, choosing the nodes — not just the weights — is precisely the freedom that orthogonality converts into $n$ extra exact degrees.

Exercises Intermediate+

Exercise 3 (medium, symbolic).

Derive the two-point Gauss-Legendre weights from the interpolatory formula $w_{i} = \int_{- 1}^{1} L_{i} d x$ on nodes $x_{1, 2} = \mp 1/ 3$ , and confirm they are equal and positive.

Hint

$L_{1} (x) = (x - x_{2}) / (x_{1} - x_{2})$ and $L_{2} (x) = (x - x_{1}) / (x_{2} - x_{1})$ , with $x_{1} - x_{2} = - 2/ 3$ . Use $\int_{- 1}^{1} 1 d x = 2$ and $\int_{- 1}^{1} x d x = 0$ .

Answer

With $x_{1} = - 1/ 3$ , $x_{2} = + 1/ 3$ : $w_{1} = \int_{- 1}^{1} (x - x_{2}) / (x_{1} - x_{2}) d x = \frac{1}{x _{1} - x _{2}} (\int_{- 1}^{1} x d x - x_{2} \int_{- 1}^{1} 1 d x) = \frac{0 - 2 x _{2}}{x _{1} - x _{2}}$ . Substituting $x_{1} - x_{2} = - 2/ 3$ and $x_{2} = 1/ 3$ : $w_{1} = \frac{- 2/ 3}{- 2/ 3} = 1$ . By the symmetric computation $w_{2} = 1$ . Both are positive and sum to $2 = \int_{- 1}^{1} 1 d x$ , confirming the corollary. Rubric: full credit for both cardinal integrals and the equality $w_{1} = w_{2} = 1$ .

Exercise 8 (hard, symbolic).

Derive the Gauss error term $E (f) = \frac{( ϕ _{n} , ϕ _{n} ) _{w}}{k _{n}^{2} ( 2 n )!} f^{(2 n)} (η)$ for $f \in C^{2 n} [a, b]$ , where $k_{n}$ is the leading coefficient of $ϕ_{n}$ . Use Hermite interpolation at the doubled nodes.

Hint

Let $H \in Π_{2 n - 1}$ be the Hermite interpolant matching $f$ and $f^{'}$ at the $n$ nodes. The rule integrates $H$ exactly, and $f - H = \frac{f ^{(2 n)} ( ξ )}{( 2 n )!} \prod_{i} (x - x_{i})^{2}$ ; note $\prod_{i} (x - x_{i}) = ϕ_{n} / k_{n}$ .

Answer

Let $H \in Π_{2 n - 1}$ be the Hermite interpolant of $f$ at the doubled nodes, so $H (x_{i}) = f (x_{i})$ and $H^{'} (x_{i}) = f^{'} (x_{i})$ . The Hermite error formula gives $f (x) - H (x) = \frac{f ^{(2 n)} ( ξ ( x ))}{( 2 n )!} \prod_{i = 1}^{n} (x - x_{i})^{2}$ for some $ξ (x)$ . Since the Gauss rule is exact on $Π_{2 n - 1}$ , $Q (f) = Q (H) = \int_{a}^{b} w H d x$ (and $Q$ uses only values $f (x_{i}) = H (x_{i})$ ). Therefore $$ E(f) = \int_a^b w f - Q(f) = \int_a^b w(f - H),dx = \frac{1}{(2n)!}\int_a^b w(x),f^{(2n)}(\xi(x))\prod_i (x-x_i)^2,dx . $$ The factor $\prod_{i} (x - x_{i})^{2} = (ϕ_{n} / k_{n})^{2} \geq 0$ is a fixed nonnegative weight, so the weighted mean value theorem supplies $η$ with $E (f) = \frac{f ^{(2 n)} ( η )}{( 2 n )!} \int_{a}^{b} w (ϕ_{n} / k_{n})^{2} d x = \frac{f ^{(2 n)} ( η )}{( 2 n )! k _{n}^{2}} (ϕ_{n}, ϕ_{n})_{w}$ . Rubric: full credit for the Hermite interpolant, the exactness step, the $(ϕ_{n} / k_{n})^{2}$ identification, and the mean-value collapse.

Advanced results Master

The elementary theory fixes the nodes and weights and proves the exactness ceiling; the master layer turns the construction into a spectral algorithm, controls the error and convergence precisely, and places the Radau and Lobatto variants inside the same orthogonal-polynomial frame.

Theorem 1 (the error term and its sign). For $f \in C^{2 n} [a, b]$ the Gauss rule has error $$ E(f) = \frac{(\phi_n, \phi_n)_w}{k_n^2,(2n)!},f^{(2n)}(\eta), \qquad \eta \in (a,b), $$ where $k_{n}$ is the leading coefficient of $ϕ_{n}$ , obtained by integrating the Hermite-interpolation remainder $f - H = f^{(2 n)} (ξ) / (2 n)! \prod_{i} (x - x_{i})^{2}$ against $w$ and collapsing the one-signed factor $(ϕ_{n} / k_{n})^{2}$ by the weighted mean value theorem ^{[Süli-Mayers 2003 §10.4]}. The error is governed by the $2 n$ -th derivative — twice the order an $n$ -point interpolatory rule reaches — and the constant is the squared $L^{2}$ norm of the monic orthogonal polynomial. For Gauss-Legendre this evaluates to $E (f) = \frac{2 ^{2 n + 1} ( n ! ) ^{4}}{( 2 n + 1 ) [( 2 n )! ] ^{3}} f^{(2 n)} (η)$ , decaying super-geometrically in $n$ for analytic $f$ .

Theorem 2 (convergence for every continuous integrand). For every $f \in C [a, b]$ the Gauss rules converge, $Q_{n} (f) \to \int_{a}^{b} w f d x$ as $n \to \infty$ . The positivity corollary gives $\sum_{i} ∣ w_{i}^{(n)} ∣ = \sum_{i} w_{i}^{(n)} = \int_{a}^{b} w$ for all $n$ , a uniform bound on the absolute-weight sums; since the rules are exact on a polynomial of each degree and polynomials are dense in $(C [a, b], ∥ \cdot ∥_{\infty})$ by Weierstrass, the Pólya-Steklov criterion of 43.09.01 applies and forces convergence ^{[Davis-Rabinowitz Ch. 3]}. This is the structural reward of positivity: Gauss rules converge for all continuous functions, where the equispaced Newton-Cotes family diverges, and the divergence-versus-convergence dichotomy is decided entirely by the boundedness of $\sum ∣ w_{i}^{(n)} ∣$ .

Theorem 3 (Golub-Welsch: nodes and weights as a symmetric eigenproblem). In the orthonormal scaling $\hat{ϕ}_{k} = ϕ_{k} /∥ ϕ_{k} ∥_{w}$ , the three-term recurrence of 43.08.05 reads $x \hat{ϕ}_{k} = β_{k + 1} \hat{ϕ}_{k + 1} + α_{k} \hat{ϕ}_{k} + β_{k} \hat{ϕ}_{k - 1}$ , so multiplication by $x$ on $span {\hat{ϕ}_{0}, \dots, \hat{ϕ}_{n - 1}}$ is the symmetric tridiagonal Jacobi matrix $$ J_n = \begin{pmatrix} \alpha_0 & \sqrt{\beta_1} & & \ \sqrt{\beta_1} & \alpha_1 & \sqrt{\beta_2} & \ & \ddots & \ddots & \ddots \ & & \sqrt{\beta_{n-1}} & \alpha_{n-1} \end{pmatrix}. $$ Its eigenvalues are exactly the Gauss nodes (the roots of $ϕ_{n}$ ), and if $v^{(i)}$ is the normalised eigenvector for node $x_{i}$ , the Gauss weight is $w_{i} = β_{0} (v_{1}^{(i)})^{2}$ with $β_{0} = \int_{a}^{b} w d x$ — the squared first component scaled by the total mass. Thus computing an $n$ -point Gauss rule reduces to one symmetric tridiagonal eigendecomposition, solved in $O (n^{2})$ by the symmetric QR/QL iteration, which is the Golub-Welsch algorithm ^{[Gautschi 2004]}. This routes the entire construction through the symmetric eigenproblem of 43.06.01: real eigenvalues and an orthonormal eigenbasis are guaranteed by symmetry, and the simplicity of the nodes is the simplicity of a tridiagonal matrix with nonzero off-diagonals.

Theorem 4 (Radau and Lobatto from modified weights). Prescribing nodes trades exactness for endpoint control, and each variant is itself a Gauss rule for a modified weight. The Gauss-Radau rule fixes one endpoint, say $a$ : its free nodes are the roots of the degree- $(n - 1)$ orthogonal polynomial for the modified weight $(x - a) w (x)$ , and it is exact on $Π_{2 n - 2}$ . The Gauss-Lobatto rule fixes both endpoints: its $n - 2$ free nodes are the roots of the orthogonal polynomial for $(x - a) (b - x) w (x)$ , and it is exact on $Π_{2 n - 3}$ ^{[Gautschi 2004]}. The weights remain positive by the same $\int w L_{i}^{2}$ argument applied to the reduced free-node set together with explicit positive endpoint weights, and at the matrix level Radau and Lobatto are rank-one and rank-two modifications of $J_{n}$ that plant the prescribed eigenvalues. Lobatto rules underlie the spectral-element and collocation methods that need function values at element boundaries.

Synthesis. The division identity $p = q ϕ_{n} + r$ is the foundational reason Gauss quadrature, orthogonal polynomials, and the symmetric eigenproblem are one subject: exactness on $Π_{2 n - 1}$ is equivalent to the node polynomial being orthogonal to $Π_{n - 1}$ , which forces the nodes to be the roots of $ϕ_{n}$ , which are the eigenvalues of the Jacobi matrix. The central insight is that the freedom to choose nodes — the freedom the equispaced rules discard — is exactly the orthogonality that doubles the exactness, and that the same orthogonality forces the weights positive, so the very property that maximises accuracy also guarantees stability and convergence; this is exactly the dichotomy of 43.09.01, where sign-definiteness and bounded absolute-weight sums separate the usable rules from the divergent ones, now resolved in Gauss's favour by construction. Putting these together, the Gauss rule is dual to the least-squares projection of 43.08.05: there $ϕ_{n}$ is the residual direction perpendicular to $Π_{n - 1}$ under the integral, here its roots are the optimal sampling points, and the bridge between them is the single Gram-Schmidt-against-a-weight construction that simultaneously diagonalises the least-squares normal equations, generates the recurrence coefficients of the Jacobi matrix, and hands the Golub-Welsch algorithm its eigenproblem. This generalises the interpolatory quadrature of 43.09.01 from degree $n - 1$ to $2 n - 1$ and is the integration-theoretic shadow of the spectral theorem of 43.06.01, one tridiagonal matrix governing both where to sample and how much each sample counts.

Full proof set Master

Proposition 1 (maximal exactness of the Gauss rule). The interpolatory rule on the roots of $ϕ_{n}$ is exact on $Π_{2 n - 1}$ , and no $n$ -point rule is exact on $Π_{2 n}$ .

Proof. For $p \in Π_{2 n - 1}$ , polynomial division by $ϕ_{n}$ (degree $n$ ) gives unique $q, r \in Π_{n - 1}$ with $p = q ϕ_{n} + r$ . Then $\int_{a}^{b} w p = (ϕ_{n}, q)_{w} + \int_{a}^{b} w r = \int_{a}^{b} w r$ , since $q \in Π_{n - 1}$ and $ϕ_{n} ⊥ Π_{n - 1}$ by 43.08.05. At each node $x_{i}$ , $ϕ_{n} (x_{i}) = 0$ , so $p (x_{i}) = r (x_{i})$ and $Q (p) = Q (r)$ ; the interpolatory rule on $n$ nodes is exact on $Π_{n - 1} ∋ r$ , so $Q (r) = \int_{a}^{b} w r$ . Hence $Q (p) = \int_{a}^{b} w p$ . For the ceiling, $g = \prod_{i} (x - x_{i})^{2} \in Π_{2 n}$ has $Q (g) = 0$ (vanishes at all nodes) while $\int_{a}^{b} w g > 0$ ( $g > 0$ off the finite node set, $w > 0$ ), so $E (g) \neq = 0$ for any $n$ -point rule. $□$

Proposition 2 (the orthogonal-polynomial nodes are forced). If an $n$ -point rule is exact on $Π_{2 n - 1}$ , its nodes are the roots of $ϕ_{n}$ .

Proof. Let $y_{1}, \dots, y_{n}$ be the nodes and $ψ (x) = \prod_{i} (x - y_{i}) \in Π_{n}$ , monic of degree $n$ . For arbitrary $s \in Π_{n - 1}$ , $s ψ \in Π_{2 n - 1}$ vanishes at every node, so exactness gives $0 = Q (s ψ) = \int_{a}^{b} w s ψ d x = (ψ, s)_{w}$ . Thus $ψ ⊥ Π_{n - 1}$ . A monic degree- $n$ polynomial orthogonal to $Π_{n - 1}$ equals the monic $ϕ_{n}$ : their difference lies in $Π_{n - 1}$ and is orthogonal to $Π_{n - 1}$ , hence to itself, hence zero. So $ψ = ϕ_{n}$ (monic scaling) and the $y_{i}$ are its roots. $□$

Proposition 3 (positivity and total mass of the weights). The Gauss weights satisfy $w_{i} = \int_{a}^{b} w L_{i}^{2} d x > 0$ and $\sum_{i} w_{i} = \int_{a}^{b} w d x$ .

Proof. $L_{i}^{2} \in Π_{2 n - 2} \subset Π_{2 n - 1}$ , so exactness gives $\int_{a}^{b} w L_{i}^{2} = Q (L_{i}^{2}) = \sum_{j} w_{j} L_{i} (x_{j})^{2} = w_{i}$ , using $L_{i} (x_{j}) = δ_{ij}$ . Since $w > 0$ and $L_{i}$ is a nonzero polynomial, $\int_{a}^{b} w L_{i}^{2} > 0$ , so $w_{i} > 0$ . Applying exactness to $1 \in Π_{0}$ gives $\sum_{i} w_{i} = \int_{a}^{b} w d x$ . $□$

Proposition 4 (the error term). For $f \in C^{2 n} [a, b]$ , $E (f) = \frac{( ϕ _{n} , ϕ _{n} ) _{w}}{k _{n}^{2} ( 2 n )!} f^{(2 n)} (η)$ for some $η \in (a, b)$ , with $k_{n}$ the leading coefficient of $ϕ_{n}$ .

Proof. Let $H \in Π_{2 n - 1}$ be the Hermite interpolant of $f$ at the doubled nodes, matching $f$ and $f^{'}$ at each $x_{i}$ . The Hermite remainder is $f (x) - H (x) = \frac{f ^{(2 n)} ( ξ ( x ))}{( 2 n )!} ω (x)^{2}$ with $ω (x) = \prod_{i} (x - x_{i}) = ϕ_{n} (x) / k_{n}$ . Exactness on $Π_{2 n - 1}$ gives $Q (f) = Q (H) = \int_{a}^{b} w H$ , since $Q$ depends only on the values $f (x_{i}) = H (x_{i})$ . Hence $E (f) = \int_{a}^{b} w (f - H) = \frac{1}{( 2 n )!} \int_{a}^{b} w f^{(2 n)} (ξ (x)) ω (x)^{2} d x$ . The factor $w ω^{2} \geq 0$ and $f^{(2 n)}$ is continuous, so the weighted mean value theorem for integrals yields $η$ with $E (f) = \frac{f ^{(2 n)} ( η )}{( 2 n )!} \int_{a}^{b} w ω^{2} = \frac{f ^{(2 n)} ( η )}{( 2 n )! k _{n}^{2}} \int_{a}^{b} w ϕ_{n}^{2} = \frac{( ϕ _{n} , ϕ _{n} ) _{w}}{k _{n}^{2} ( 2 n )!} f^{(2 n)} (η)$ . $□$

Proposition 5 (Golub-Welsch eigenvalue characterisation). The Gauss nodes are the eigenvalues of the symmetric tridiagonal Jacobi matrix $J_{n}$ , and $w_{i} = β_{0} (v_{1}^{(i)})^{2}$ where $v^{(i)}$ is the normalised eigenvector at node $x_{i}$ and $β_{0} = \int_{a}^{b} w$ .

Proof. Write $\hat{ϕ} (x) = (\hat{ϕ}_{0} (x), \dots, \hat{ϕ}_{n - 1} (x))^{⊤}$ . The orthonormal recurrence gives $x \hat{ϕ} (x) = J_{n} \hat{ϕ} (x) + β_{n} \hat{ϕ}_{n} (x) e_{n}$ . At a root $x_{i}$ of $ϕ_{n}$ (equivalently of $\hat{ϕ}_{n}$ ), the trailing term drops and $J_{n} \hat{ϕ} (x_{i}) = x_{i} \hat{ϕ} (x_{i})$ , so $x_{i}$ is an eigenvalue of $J_{n}$ with eigenvector $\hat{ϕ} (x_{i})$ (nonzero since $\hat{ϕ}_{0} \equiv 1/ β_{0} \neq = 0$ ). The $n$ distinct roots give $n$ distinct eigenvalues, exhausting the spectrum of the $n \times n$ symmetric $J_{n}$ , whose eigenvalues are real and whose eigenvectors are orthogonal by the spectral theorem of 43.06.01. The normalised eigenvector is $v^{(i)} = \hat{ϕ} (x_{i}) /∥ \hat{ϕ} (x_{i}) ∥$ , with first component $v_{1}^{(i)} = \hat{ϕ}_{0} (x_{i}) /∥ \hat{ϕ} (x_{i}) ∥ = (1/ β_{0}) /∥ \hat{ϕ} (x_{i}) ∥$ . The Christoffel-Darboux reproducing-kernel weight formula $w_{i} = 1/ \sum_{k = 0}^{n - 1} \hat{ϕ}_{k} (x_{i})^{2} = 1/∥ \hat{ϕ} (x_{i}) ∥^{2}$ of 43.08.05 then gives $β_{0} (v_{1}^{(i)})^{2} = β_{0} \cdot \frac{1/ β _{0}}{∥ ϕ ^ ( x _{i} ) ∥ ^{2}} = \frac{1}{∥ ϕ ^ ( x _{i} ) ∥ ^{2}} = w_{i}$ . $□$

Proposition 6 (Radau exactness from the modified weight). The Gauss-Radau rule with fixed node $a$ — free nodes the roots of the degree- $(n - 1)$ orthogonal polynomial for $(x - a) w (x)$ — is exact on $Π_{2 n - 2}$ .

Proof. Let $\tilde{w} (x) = (x - a) w (x) \geq 0$ on $[a, b]$ , and let $\tilde{ϕ}_{n - 1}$ be its degree- $(n - 1)$ orthogonal polynomial with roots $x_{2}, \dots, x_{n} \in (a, b)$ ; the Radau nodes are $x_{1} = a$ together with these. Take $p \in Π_{2 n - 2}$ and divide by $(x - a) \tilde{ϕ}_{n - 1} (x)$ , degree $n$ : $p = q \cdot (x - a) \tilde{ϕ}_{n - 1} + r$ with $de g q \leq n - 2$ , $de g r \leq n - 1$ . Integrating against $w$ : $\int_{a}^{b} w p = \int_{a}^{b} w q (x - a) \tilde{ϕ}_{n - 1} + \int_{a}^{b} w r = \int_{a}^{b} \tilde{w} q \tilde{ϕ}_{n - 1} + \int_{a}^{b} w r = (\tilde{ϕ}_{n - 1}, q)_{\tilde{w}} + \int_{a}^{b} w r$ . Since $q \in Π_{n - 2}$ and $\tilde{ϕ}_{n - 1} ⊥ Π_{n - 2}$ under $(\cdot, \cdot)_{\tilde{w}}$ , the first term vanishes, leaving $\int_{a}^{b} w p = \int_{a}^{b} w r$ . The rule has $n$ nodes and interpolatory weights, exact on $Π_{n - 1} ∋ r$ ; and at each node the factor $(x - a) \tilde{ϕ}_{n - 1}$ vanishes (at $a$ by $(x - a)$ , at the others by $\tilde{ϕ}_{n - 1}$ ), so $p (x_{i}) = r (x_{i})$ and $Q (p) = Q (r) = \int_{a}^{b} w r = \int_{a}^{b} w p$ . $□$

Connections Master

The Gauss nodes are the roots of the orthogonal polynomial $ϕ_{n}$ constructed in 43.08.05 — real, simple, interlacing, and inside $(a, b)$ — and the weight-positivity, error constant $(ϕ_{n}, ϕ_{n})_{w} / k_{n}^{2}$ , and Golub-Welsch eigenvector formula all descend from the three-term recurrence and Christoffel-Darboux kernel developed there; this unit is the quadrature payoff of that unit's root-structure theorem, the forward reference its own Synthesis anticipated.
The maximal-exactness theorem is the resolution of the degree-of-exactness ceiling exposed in 43.09.01: where the interpolatory Newton-Cotes rules on $n$ fixed equispaced nodes reach only degree $n - 1$ (or $n$ by symmetry) and acquire fatal negative weights at high order, the Gauss rule on the orthogonal-polynomial nodes reaches $2 n - 1$ with strictly positive weights summing to $\int w$ , so the Pólya-Steklov convergence criterion proved there now certifies convergence for every continuous integrand rather than divergence.
The Golub-Welsch algorithm routes node-and-weight computation through the symmetric tridiagonal eigenproblem of 43.06.01: the Jacobi matrix of recurrence coefficients is real symmetric, so its eigenvalues (the nodes) are real and its eigenvectors orthogonal by the spectral theorem, and the symmetric QR/QL iteration computes the full rule in $O (n^{2})$ ; the simplicity of the Gauss nodes is the simplicity of a tridiagonal matrix with nonzero subdiagonal.
Gauss quadrature is the integration tool behind spectral and pseudospectral discretisations and the Gaussian-process / random-matrix computations of 37.08.01, where weighted integrals against $e^{- x^{2}}$ (Gauss-Hermite) and against equilibrium measures are evaluated to machine precision; the Lobatto variant of Theorem 4, fixing both endpoints, supplies the boundary collocation nodes used in the spectral-element methods that border the finite-element theory of 24.04.07.
The error term's reliance on the Hermite interpolant at doubled nodes ties this unit to the Hermite and osculatory interpolation of 43.08.03, and the weighted mean value theorem that collapses the one-signed factor $(ϕ_{n} / k_{n})^{2}$ is the same integral mean value theorem underlying the Lagrange remainder of 02.05.02, so the Gauss error is the $2 n$ -th-order member of the single family of derivative-times-fixed-kernel error formulas that also produced the Peano representation of 43.09.01.

Historical & philosophical context Master

Carl Friedrich Gauss introduced the method in Methodus nova integralium valores per approximationem inveniendi, read to the Göttingen Society in 1814 and published in its Commentationes the same year. Gauss derived the $n$ -point rule on $[- 1, 1]$ exact to degree $2 n - 1$ through the continued-fraction expansion of a logarithmic generating function — equivalently, through what are now recognised as the Legendre polynomials — and tabulated nodes and weights for $n$ up to $7$ , doubling the attainable exactness of the equispaced interpolatory rules of Newton and Cotes for the same number of evaluations ^{[Gauss 1814]}.

The orthogonal-polynomial reading that this unit follows is due to Carl Gustav Jacobi, who in an 1826 paper recast Gauss's nodes as the roots of the Legendre polynomials and supplied the modern proof via the division argument, freeing the construction from continued fractions. The generalisation to an arbitrary positive weight is the work of Pafnuty Chebyshev and Elwin Bruno Christoffel — the weights $w_{i}$ carry Christoffel's name as the Christoffel numbers from his 1858 study — and of Thomas Joannes Stieltjes, whose 1884 work tied the weights and nodes to the moment problem and continued fractions. The Radau and Lobatto endpoint variants are due to Rodolphe Radau (1880) and Rehuel Lobatto (1852). The reduction of the entire computation to a symmetric tridiagonal eigenproblem was published by Gene Golub and John Welsch in 1969, and the recurrence-coefficient viewpoint that makes it practical was developed by Walter Gautschi ^{[Gautschi 1981]}.

Bibliography Master

@book{sulimayers2003,
  author    = {S\"{u}li, Endre and Mayers, David F.},
  title     = {An Introduction to Numerical Analysis},
  publisher = {Cambridge University Press},
  year      = {2003}
}

@book{davisrabinowitz1984,
  author    = {Davis, Philip J. and Rabinowitz, Philip},
  title     = {Methods of Numerical Integration},
  edition   = {2},
  publisher = {Academic Press},
  year      = {1984}
}

@book{gautschi2004,
  author    = {Gautschi, Walter},
  title     = {Orthogonal Polynomials: Computation and Approximation},
  publisher = {Oxford University Press},
  year      = {2004}
}

@article{golubwelsch1969,
  author  = {Golub, Gene H. and Welsch, John H.},
  title   = {Calculation of {G}auss Quadrature Rules},
  journal = {Mathematics of Computation},
  volume  = {23},
  number  = {106},
  year    = {1969},
  pages   = {221--230}
}

@article{gauss1814,
  author  = {Gauss, Carl Friedrich},
  title   = {Methodus nova integralium valores per approximationem inveniendi},
  journal = {Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores},
  volume  = {3},
  year    = {1814},
  pages   = {39--76}
}

@incollection{gautschi1981,
  author    = {Gautschi, Walter},
  title     = {A Survey of {G}auss-{C}hristoffel Quadrature Formulae},
  booktitle = {E. B. Christoffel: The Influence of His Work on Mathematics and the Physical Sciences},
  editor    = {Butzer, P. L. and Feh\'{e}r, F.},
  publisher = {Birkh\"{a}user},
  year      = {1981},
  pages     = {72--147}
}

@article{jacobi1826,
  author  = {Jacobi, Carl Gustav Jacob},
  title   = {Ueber {G}au{\ss}' neue Methode, die Werthe der Integrale n\"{a}herungsweise zu finden},
  journal = {Journal f\"{u}r die reine und angewandte Mathematik},
  volume  = {1},
  year    = {1826},
  pages   = {301--308}
}

Prerequisites

43.09.01
43.08.05

Tier anchors

beginner: Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §10.1-10.3 (Gauss quadrature as the choice of both nodes and weights to maximise exactness, at a first-course level, with the two-point Gauss-Legendre rule worked numerically); Burden-Faires 2011 *Numerical Analysis* 9e (Brooks/Cole) §4.7 (Gaussian quadrature, the Legendre nodes and weights tabulated and applied to worked integrals)
intermediate: Süli-Mayers 2003 *An Introduction to Numerical Analysis* (Cambridge) §10.1-10.4 (the 2n-1 exactness theorem with nodes at the orthogonal-polynomial roots, positivity of the weights, and the f^{(2n)} error term); Stoer-Bulirsch 2002 *Introduction to Numerical Analysis* 3e (Springer) §3.6 (Gaussian integration, the orthogonal-polynomial nodes, the weight formula, and the error representation)
master: Davis-Rabinowitz 1984 *Methods of Numerical Integration* 2e (Academic Press) Ch. 2-3 (Gauss rules, the orthogonal-polynomial connection, the Radau and Lobatto variants, and the Stieltjes error theory); Gautschi 2004 *Orthogonal Polynomials: Computation and Approximation* (Oxford) Ch. 1-3 (the Jacobi matrix, the Golub-Welsch eigenvalue algorithm, and the classical Gauss families); Gautschi 1981 'A survey of Gauss-Christoffel quadrature formulae' in *E. B. Christoffel* (Birkhäuser) (the comprehensive survey of Gauss-Christoffel rules, the Radau/Lobatto modifications, and the weight-positivity and convergence theory)

References

Süli, E. & Mayers, D. F. — An Introduction to Numerical Analysis · Cambridge University Press (2003). Chapter 10 ('Numerical integration -- II') develops Gauss quadrature: §10.1 poses the optimal-quadrature problem of choosing the n nodes and n weights of a rule Q(f)=sum_{i=1}^n w_i f(x_i) approximating int_a^b w(x) f(x) dx to integrate exactly polynomials of the highest possible degree, and shows by a parameter count (2n free parameters) that 2n-1 is the ceiling. §10.2 proves the central theorem: the rule is exact on Pi_{2n-1} if and only if its nodes are the n roots of the degree-n orthogonal polynomial phi_n for the weight w and its weights are the interpolatory weights w_i=int_a^b w L_i; the proof divides an arbitrary p in Pi_{2n-1} as p = q phi_n + r with q,r in Pi_{n-1} and uses (phi_n,q)_w=0. §10.3 proves the weights are strictly positive (w_i = int w L_i^2 > 0, since L_i^2 in Pi_{2n-2} is integrated exactly) and derives the convergence of Gauss rules for every continuous f from positivity and the Weierstrass theorem. §10.4 gives the error term E(f)=f^{(2n)}(eta)/(2n)! times (phi_n,phi_n)_w / k_n^2 (with k_n the leading coefficient), the Gauss-Legendre, Gauss-Chebyshev, Gauss-Hermite and Gauss-Laguerre families as the classical weights, and notes the Gauss-Radau and Gauss-Lobatto variants that prescribe one or both endpoints.
Davis, P. J. & Rabinowitz, P. — Methods of Numerical Integration · Academic Press, 2nd edition (1984). Chapters 2-3 give the comprehensive account of Gauss-Christoffel quadrature: the orthogonal-polynomial characterisation of the nodes, the positivity of the Christoffel numbers (weights) w_i = int w lambda_i with the closed form w_i = -k_{n+1}/(k_n) (phi_n,phi_n)_w / (phi_n'(x_i) phi_{n+1}(x_i)) and the reproducing-kernel form w_i = 1/sum_{k=0}^{n-1} hat-phi_k(x_i)^2, the f^{(2n)} error term derived from Hermite interpolation at the doubled nodes, the Gauss-Radau rule (one fixed endpoint, exact on Pi_{2n-2}) and the Gauss-Lobatto rule (both endpoints fixed, exact on Pi_{2n-3}) obtained from orthogonal polynomials with respect to the modified weights (x-a)w, (b-x)w and (x-a)(b-x)w, and the convergence theorem for continuous integrands via uniform boundedness of the positive weights summing to int w.
Gautschi, W. — Orthogonal Polynomials: Computation and Approximation · Oxford University Press (2004). Chapters 1-3 take the recurrence coefficients (alpha_k, beta_k) of the orthogonal family as the computational object and reduce Gauss-rule construction to a symmetric tridiagonal eigenproblem: the Jacobi matrix J_n (alpha_k on the diagonal, sqrt(beta_k) on the off-diagonals) has the Gauss nodes as its eigenvalues and the Gauss weights as beta_0 = int w times the squares of the first components of the normalised eigenvectors -- the Golub-Welsch algorithm, computing an n-point Gauss rule in O(n^2) operations via the symmetric tridiagonal QR/QL iteration. The chapters also tabulate the classical recurrence coefficients (Legendre, Chebyshev, Hermite, Laguerre, Jacobi) and treat the Radau and Lobatto modifications as rank-one or rank-two changes to J_n that fix prescribed nodes.
Gauss, C. F. — Methodus nova integralium valores per approximationem inveniendi · Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores 3 (1814), 39-76. The original memoir introducing Gaussian quadrature: Gauss derives, via continued fractions of the logarithmic series (equivalently the Legendre orthogonal polynomials), the n-point rule on [-1,1] that integrates exactly all polynomials of degree up to 2n-1, tabulates the nodes and weights for n up to 7, and demonstrates the doubling of the attainable exactness over equispaced interpolatory rules with the same number of evaluations.

Estimated time

beginner: 20m
intermediate: 50m
master: 90m