Prev
Unit 8 · Numerical Differentiation and Integration
Newton–Cotes and Composite Rules
55m100 XP
Next
∫Numerical Differentiation and Integration

Newton–Cotes and Composite Rules

Trapezoid and Simpson rules, degree of precision, composite rules and their error, total impulse

Linz, autumn 1613
Johannes Kepler, Imperial Mathematician and author of the laws of planetary motion, is getting married. He buys a supply of wine for the wedding, and the merchant comes to measure how much each barrel holds. He pushes a graduated rod through the bung hole at the top of the barrel diagonally down to the far edge of the bottom, reads off a number, and charges accordingly — whatever the barrel's actual shape. Kepler is not satisfied. One diagonal cannot capture the volume of barrels that bulge differently. He sets out to find the volumes of solids of revolution, slicing them into thin discs and summing, and in 1615 publishes Nova stereometria doliorum vinariorum — a new solid geometry of wine barrels. Among its results is a rule for a barrel's volume from three cross-sections — the two ends and the middle: V≈L6(Aend+4Amiddle+Aend)V \approx \frac{L}{6}(A_{\text{end}} + 4A_{\text{middle}} + A_{\text{end}}). German textbooks still call it Keplersche Fassregel, Kepler's barrel rule. English ones call it Simpson's rule, after Thomas Simpson, who published it 128 years later.
To integrate ff over [a,b][a, b], replace it by its interpolating polynomial on equally spaced points (Unit 7) and integrate that exactly. These are the Newton–Cotes rules, after Newton and Roger Cotes, whose tables appeared in 1722. The first two:
∫abf dx≈b−a2(f(a)+f(b)),∫abf dx≈b−a6(f(a)+4f(a+b2)+f(b))\int_a^b f\,dx \approx \frac{b-a}{2}\big(f(a)+f(b)\big), \qquad \int_a^b f\,dx \approx \frac{b-a}{6}\Big(f(a)+4f\big(\tfrac{a+b}{2}\big)+f(b)\Big)

The trapezoid rule (a straight line through two points) and Simpson's rule (a parabola through three). Kepler's barrel rule is Simpson's rule with the cross-sectional area as the integrand.

Derivation
Finding the weights, and Simpson's free extra degree
A rule of the form ∫−hhf dx≈w0f(−h)+w1f(0)+w2f(h)\int_{-h}^{h}f\,dx \approx w_0f(-h) + w_1f(0) + w_2f(h) has three unknown weights. Ask it to be exact for 11, xx and x2x^2 — three equations:
1: w0+w1+w2=2h,x: −hw0+hw2=0,x2: h2w0+h2w2=23h3.1:\ w_0 + w_1 + w_2 = 2h, \qquad x:\ -hw_0 + hw_2 = 0, \qquad x^2:\ h^2w_0 + h^2w_2 = \tfrac23h^3.
The second gives w0=w2w_0 = w_2, the third w0=w2=h/3w_0 = w_2 = h/3, and the first w1=4h/3w_1 = 4h/3:
∫−hhf dx≈h3(f(−h)+4f(0)+f(h)),\int_{-h}^{h}f\,dx \approx \frac h3\big(f(-h) + 4f(0) + f(h)\big),
which, with b−a=2hb - a = 2h, is Simpson's rule. (Two points and exactness for 11 and xx give the trapezoid weights h/2,h/2h/2, h/2 the same way.)
Degree 3 for free. For x3x^3 both sides are zero — the integral of an odd function over a symmetric interval, and the rule's −h3+h3-h^3 + h^3. The first power it misses is x4x^4: the rule gives h3⋅2h4=23h5\tfrac h3 \cdot 2h^4 = \tfrac23h^5, the integral 25h5\tfrac25h^5.
The error constant from that miss. For a smooth ff the error is dominated by its x4x^4 part, f(4)x4/24f^{(4)}x^4/24. The rule overestimates ∫x4\int x^4 by 23h5−25h5=415h5\tfrac23h^5 - \tfrac25h^5 = \tfrac4{15}h^5, so the error of one panel is
∫−hhf dx−S=−f(4)24⋅415h5=−h590f(4)(ξ),\int_{-h}^{h}f\,dx - S = -\frac{f^{(4)}}{24}\cdot\frac4{15}h^5 = -\frac{h^5}{90}f^{(4)}(\xi),
exactly the formula quoted below (a careful version replaces the constant f(4)f^{(4)} by its value at some ξ\xi).
💡Degree of precision
A rule has degree of precision dd if it integrates every polynomial of degree ≤d\le d exactly. The trapezoid rule has d=1d = 1. Simpson's rule, built from a parabola, should have d=2d = 2 — but by symmetry the error of the cubic term cancels too, and it has d=3d = 3: on [0,1][0, 1] it gives exactly 1/31/3 for x2x^2 and 1/41/4 for x3x^3, but 0.20830.2083 instead of 0.20.2 for x4x^4. The error of one Simpson panel of width 2h2h is −h590f(4)(ξ)-\tfrac{h^5}{90}f^{(4)}(\xi).
High-degree Newton–Cotes rules on one panel suffer exactly Runge's problem (lesson 7-2): from nine points on, some weights become negative and errors in ff are amplified. The practical answer is composite rules: split [a,b][a, b] into nn panels of width hh and apply a low-order rule on each.
Tn=h(12f0+f1+⋯+fn−1+12fn),Sn=h3(f0+4f1+2f2+4f3+⋯+4fn−1+fn)T_n = h\Big(\tfrac12f_0 + f_1 + \dots + f_{n-1} + \tfrac12f_n\Big),\qquad S_n = \frac{h}{3}\big(f_0 + 4f_1 + 2f_2 + 4f_3 + \dots + 4f_{n-1} + f_n\big)

Composite trapezoid and Simpson rules (n even for Simpson). Errors: −(b − a)h²f″/12 and −(b − a)h⁴f⁽⁴⁾/180. Doubling n divides the trapezoid error by 4 and Simpson's by 16.

Derivation
The composite error, from one panel to many
One trapezoid panel. The trapezoid rule integrates the straight line pp through the two ends. Lesson 7-1's error formula with two nodes says f(x)−p(x)=12f′′(ξx)(x−a)(x−b)f(x) - p(x) = \tfrac12 f''(\xi_x)(x - a)(x - b). Integrate over the panel of width h=b−ah = b - a:
∫ab(x−a)(x−b) dx=−h36.\int_a^b (x - a)(x - b)\,dx = -\frac{h^3}{6}.
The factor (x−a)(x−b)(x - a)(x - b) never changes sign on the panel, so the mean value theorem for integrals lets f′′f'' come out at a single point η\eta:
∫abf dx−T=12f′′(η)(−h36)=−h312f′′(η).\int_a^b f\,dx - T = \tfrac12f''(\eta)\Big(-\frac{h^3}{6}\Big) = -\frac{h^3}{12}f''(\eta).
Many panels. The composite rule makes this error in each of its nn panels. Adding them, and replacing the sum of nn values of f′′f'' by nn times their average (itself a value of f′′f'' somewhere):
ET=−h312∑i=1nf′′(ηi)=−h312 n f′′(η)=−(b−a) h212f′′(η),E_T = -\frac{h^3}{12}\sum_{i=1}^n f''(\eta_i) = -\frac{h^3}{12}\,n\,f''(\eta) = -\frac{(b - a)\,h^2}{12}f''(\eta),
since nh=b−anh = b - a. Each panel's error is O(h3)O(h^3), but there are n∝1/hn \propto 1/h of them: the total is O(h2)O(h^2) — the same 'one order lost to accumulation' as Euler's method in lesson 9-1. Simpson's panels of width 2h2h have error h5f(4)/90h^5f^{(4)}/90 and there are (b−a)/2h(b - a)/2h of them: (b−a)h4f(4)/180(b - a)h^4f^{(4)}/180.
📐The orders, measured
For ∫0πsin⁡x dx=2\int_0^\pi \sin x\,dx = 2: with n=4,8,16,32,64n = 4, 8, 16, 32, 64 panels the trapezoid errors are 1.0×10−11.0 \times 10^{-1}, 2.6×10−22.6 \times 10^{-2}, 6.4×10−36.4 \times 10^{-3}, 1.6×10−31.6 \times 10^{-3}, 4.0×10−44.0 \times 10^{-4} — ratios 4.03, 4.01, 4.00, 4.00. Simpson: 4.6×10−34.6 \times 10^{-3} down to 6.5×10−86.5 \times 10^{-8}, ratios 16.9, 16.2, 16.06, 16.01. Watching the ratio settle at 2p2^p is the standard check that a quadrature code — or any discretisation — is implemented correctly.
Nature uses as little as possible of anything.
– Attributed to Johannes Kepler, as quoted in the Viking Book of Aphorisms
📐Total impulse from a static fire
Thrust sampled at 1 kHz for 3 s, with a little sensor noise. The trapezoid rule gives a total impulse of 6169.0 N s; Simpson gives the same to the reported precision. With 2.90 kg of propellant, the delivered specific impulse is I/(mpg0)=216.9I/(m_pg_0) = 216.9 s. The noise changes the impulse by only 0.32 N s: integration averages random errors, the opposite of differentiation. A logger recording only 20 samples per second would have reported 6143.8 N s — off by 0.4%, because the sharp ignition rise and burnout are under-resolved. For data, the grid is given, and the trapezoid rule is usually the honest choice.
⚠Kinks spoil the order
The error formulas assume f′′f'' (or f(4)f^{(4)}) exists across each panel. A thrust curve has a kink at burnout; a stress–strain curve at yield; a piecewise material property at its breakpoints. A panel straddling the kink drops the rule to first or second order. Place panel boundaries at known kinks, or split the integral there.
💡Where you will use this
Rocket Propulsion: total and specific impulse from test-stand data. Structures: section properties (area, centroid, second moment) of shapes defined by measured coordinates; work and energy from load–displacement curves. Hydraulics: discharge from velocity profiles measured across a river (the velocity–area method). Orbital Mechanics: averaging perturbations over an orbit.
📚Further reading
  • Philip Davis and Philip Rabinowitz, Methods of Numerical Integration (2nd ed., Academic Press, 1984; Dover reprint 2007) — the encyclopedia of quadrature rules and their errors.
  • Endre Süli and David Mayers, An Introduction to Numerical Analysis (Cambridge University Press, 2003), chapter 7, 'Numerical integration — I' — Newton–Cotes and composite rules with full error proofs.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 8, 'Numerical Integration and Differentiation'.
↓But there is a problem
Composite rules spend their points evenly. For a smooth integrand that is fine. For one with a sharp feature — a peak, a near-singularity, the square-root behaviour near a crack tip — most of the points land where the function is easy, and the error is set by the few that land where it is hard. When you call scipy.integrate.quad, it does not ask you for a grid. It hands the problem to a Fortran library written in the early 1980s by four numerical analysts in Belgium, Austria and the United States — QUADPACK — which decides by itself where to look.