Gaussian quadrature

The formula

The n-point rule on the standard interval:

∫₋₁¹ f(x) dx ≈ w₁·f(x₁) + w₂·f(x₂) + … + wₙ·f(xₙ)

Reading: evaluate f at n prescribed points, multiply each value by its weight, add. No curve fitting at run time — the points and weights were solved once, centuries ago, and live in a table.

The points xᵢ are the roots of the degree-n Legendre polynomial Pₙ(x), and each weight is

wᵢ = 2 / ( (1 − xᵢ²) · [P'ₙ(xᵢ)]² )

Reading: the closer a node sits to the ends of the interval, the less it counts. All weights come out positive, which is part of why the method is numerically well behaved.

The two rules you will actually see inside FEA codes:

2-point:  x = ±1/√3 ≈ ±0.577350,           w = 1, 1          (exact through cubics)
3-point:  x = 0, ±√(3/5) ≈ ±0.774597,      w = 8/9, 5/9, 5/9 (exact through quintics)

Any other interval [a, b] maps onto the standard one with a stretch and a shift:

∫ₐᵇ f(x) dx ≈ (b−a)/2 · Σ wᵢ · f( (b−a)/2 · xᵢ + (a+b)/2 )

The error term is proportional to f⁽²ⁿ⁾(ξ) — the 2n-th derivative of the integrand somewhere in the interval. That one factor is the whole contract: if the function is smooth, the error collapses fast; if it has a kink, all bets are off.

Where you meet it

The stress query in your FEA post-processor. Every element stiffness matrix in a commercial solver is an integral over the element's volume, and the solver evaluates it by Gaussian quadrature — a 2×2 grid of points inside a linear quad, 2×2×2 inside a linear brick. Those interior locations are the Gauss points, and they are where the code actually computes strain and stress. When the post-processor offers "element centroid," "integration point," or "nodal (extrapolated)" output, it is telling you how far the number traveled from where it was really calculated.

The review board asking why two stress numbers disagree. The contour plot says 41 ksi at a node; the element table says 36 ksi. Neither is a typo. Stress at a node is extrapolated outward from the Gauss points and then usually averaged across neighboring elements; the raw computed values live inside the element. Knowing that resolves the argument in one sentence.

Every call to scipy.integrate.quad or MATLAB's integral. The adaptive engines underneath (QUADPACK lineage) run Gauss–Kronrod pairs — a 7-point Gauss rule and a 15-point extension sharing the same samples — and use the difference between the two answers as the error estimate. Your data-reduction script is running Gauss points whether you asked for them or not.

Hand-rolling an integral in a spreadsheet. Radiated power from a pattern cut, total impulse from a thrust model, a section property with a curved boundary. Five Gauss points of a smooth function routinely beat a fifty-row trapezoid column, and the five nodes fit in one cell range.

How it works

Newton–Cotes rules — trapezoid, Simpson — fix the sample points on an even grid and only get to choose the weights. That is n free parameters, enough to integrate polynomials of degree about n−1. Gauss's move was to free the sample locations too: 2n parameters, enough to nail every polynomial through degree 2n−1. The orthogonality of the Legendre polynomials is what makes the bookkeeping work — put the nodes at the roots of Pₙ and the error terms up to degree 2n−1 cancel identically.

The payoff is not subtle. Integrating eˣ over [−1, 1], a 3-point Gauss rule lands within 7·10⁻⁵ of the true answer; a trapezoid rule needs roughly a hundred points to match it. A 2-point Gauss rule — two function evaluations — beats 3-point Simpson on the same integrand. For smooth functions the error drops roughly geometrically as you add points, not polynomially.

Smoothness is the fine print. The error rides on the 2n-th derivative, so a function that is not smooth breaks the deal. Integrate √x on [0, 1]: 5 points get you three digits, and 100 points still leave an error near 10⁻⁷ — slow, grinding algebraic convergence instead of the usual free fall, because every derivative of √x blows up at zero. The fix is never "more points." Split the interval at the kink, substitute the singularity away, or switch to the Gauss variant whose weight function absorbs it: Gauss–Jacobi for endpoint singularities, Gauss–Laguerre for [0, ∞), Gauss–Hermite for the whole real line.

Two more properties worth knowing. First, all the nodes are strictly inside the interval — the rule never touches the endpoints, which is occasionally exactly what you want. Second, Gauss nodes don't nest: the 10-point rule reuses none of the 5-point rule's samples, so you can't cheaply check your answer by doubling n. That is the specific problem Kronrod's extension solves, and why library integrators use Gauss–Kronrod pairs rather than plain Gauss.

The FEA-specific gotchas cluster around how many points you use. "Full integration" means enough Gauss points to integrate the element's stiffness exactly for an undistorted shape; "reduced integration" deliberately uses fewer. Reduced integration is cheaper and often fixes overly-stiff behavior like shear locking, but a single-point quad or brick can no longer see certain bending-shaped deformations at its lone Gauss point — the notorious hourglass modes, which cost no strain energy and show up as a zig-zag mesh collapsing through itself. Every commercial code carries hourglass stabilization for exactly this reason, and every experienced analyst has a story about the run where it wasn't enough.

The mistake people make: treating extrapolated, averaged nodal stress as the ground truth and the Gauss-point values as an implementation detail. It is precisely backwards. Barlow showed in 1976 that for standard elements the Gauss points are the locations where computed stress is most accurate — they are superconvergent, better than the element's average pedigree suggests ^[5]^[6]. Nodal values are those good numbers pushed outward through a fitting step. When a peak stress at a re-entrant corner drives a margin, check what the integration-point values say before you write the number down.

History

Carl Friedrich Gauss built the method in 1814, deriving the node locations through continued-fraction manipulations of hypergeometric series — a characteristically heavy piece of machinery for what the result turned out to be ^[1]^[3]. Twelve years later, in 1826, Carl Gustav Jacobi found the clean explanation: the mysterious optimal nodes are simply the roots of the Legendre polynomials, and orthogonality does all the work ^[1]^[3]. Christoffel then widened the theory to general weight functions, which is why the general construction carries the name Gauss–Christoffel quadrature and the weights are still called Christoffel coefficients ^[2]^[4].

The computing era added two chapters. In 1964 the Soviet mathematician Alexander Kronrod solved the nesting problem — adding n+1 new points to an n-point Gauss rule so that the old samples get reused and the pair yields both an answer and an error estimate — the scheme inside QUADPACK and nearly every modern adaptive integrator ^[7]^[8]. In 1969 Gene Golub and John Welsch showed that the nodes and weights fall out of a symmetric tridiagonal eigenvalue problem, turning table lookup into a routine computation for any n and any weight function ^[1]^[3]. And in 1976 John Barlow, working on finite element stress recovery, proved the fact every stress analyst now leans on: the points where FEA stresses come out best are the Gauss points themselves ^[5]^[6].

Related tools

/tools/beam-deflection
/tools/stress-strain
/tools/section-modulus-rect
/tools/moment-of-inertia-shapes
/tools/vibration-natural-freq

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →