Numerical integration

The Earth cannot be put on scales. But from how earthquake waves pass through it, seismologists know its density at every depth: in 1981 Adam Dziewoński and Don Anderson combined these measurements into the PREM model, ten layers from the iron core to the ocean. The mass is a sum of thin spherical shells: the integral of 4πr²ρ(r) from the centre to the surface. The density in the model is given by polynomials, so we know the exact answer and can honestly measure the error of each way of computing the integral as a sum. Choose a method, move the number of evaluations and switch on “every layer on its own”.

Start with Simpson and nodes spaced evenly from the centre to the surface. With 4097 evaluations of the density it errs by 1.2·10⁻⁴, almost like the rectangles at 1.1·10⁻⁴, and as the number of nodes grows the error falls slowly and unevenly. Now switch on “every layer on its own”: the same rule on 1024 pieces errs by 5.6·10⁻¹³, and three Gauss points per layer, 30 evaluations, give all sixteen digits. The rules are the same, the data is the same; the difference is whether the method knows where the function jumps. The whole article is about this.

The idea

The interval is cut into pieces of length h, and on each the function is replaced by something that integrates exactly: a constant, a straight segment, a parabola. The areas of the pieces are added. If the function is smooth, on a small piece it almost coincides with its replacement, and the error of a piece is proportional to a high power of h. There are 1/h pieces, so the total error behaves as h^p. The number p is the order of the method: the larger it is, the faster accuracy grows with the number of evaluations.

Why not take the antiderivative? Because often there is none. The Gaussian curve e^−x²/2, 1/ln x, the length of an ellipse have no antiderivatives among elementary functions, and a function given by a table of measurements does not even have a formula. The geometry of the integral and the formulas of Newton–Leibniz and Simpson are worked through in the chapter on the integral of “Mathematics, the Queen of the Sciences”; here, what follows from it as an algorithm: what accuracy costs, how to know the error without knowing the answer, and where the formulas lie.

Three rules and their order

Let the interval [a, b] be divided into n pieces of length h = (b − a)/n, and let f_k be the value of the function at the point a + kh.

• Rectangles (the midpoint rule): on each piece the height is the value at its middle. The sum h·(f(a + h/2) + f(a + 3h/2) + …).
• Trapezoids: neighbouring points are joined by segments. The sum h·(f₀/2 + f₁ + … + f_n−1 + f_n/2).
• Simpson: a parabola is drawn through every three neighbouring points; n is even, with weights 1, 4, 2, 4, …, 2, 4, 1 times h/3.

Take an integral not expressible in elementary functions: the share of values of the normal distribution within one σ, the integral from −1 to 1 of the density e^−x²/2/√(2π). The answer 0.6826894921370859 is known from tables of the error function. The program computes it by three rules and prints the errors:

   n    rectangles   trapezoids      Simpson   Runge for Simpson
   4    5.116e-03   -1.017e-02    3.686e-04     6.786e-04
   8    1.265e-03   -2.526e-03    2.148e-05     2.314e-05
  16    3.154e-04   -6.305e-04    1.320e-06     1.344e-06
  32    7.878e-05   -1.576e-04    8.216e-08     8.253e-08
  64    1.969e-05   -3.938e-05    5.130e-09     5.136e-09

Double the number of pieces, and the error of rectangles and trapezoids falls fourfold, and of Simpson sixteenfold. These are the orders 2, 2 and 4. Notice one more thing: the trapezoids err twice as much as the rectangles and in the other direction. So two thirds of the rectangles plus one third of the trapezoids cancel the main error, and the result is exactly Simpson’s rule with half the step. Where h² comes from is seen in the exact error formula.

The theorem says more than “the error is of order h²”. It names exactly what is paid with: the second derivative. Functions with a bounded f″ are computed by trapezoids with an error ~ h², and where the derivative blows up or the function jumps, the guarantee disappears. For rectangles the same derivation gives an error of +(b − a)h²f″/24, half as large and of the opposite sign, which is what we saw in the table. For Simpson the argument is longer (integration by parts has to be done four times), and the answer is −(b − a)h⁴f⁗(η)/180: it is paid with the fourth derivative, and the rule is exact for any cubic polynomial.

On the chart at the bottom of the first scene logarithms are plotted on both axes, and the power law “error ≈ C/n^p” turns into a line with slope −p. When layers are computed separately, the lines of rectangles and trapezoids go down with slope −2, Simpson’s with slope −4, and Gauss lies at the bottom at once, at the level of rounding errors. When nodes are placed blindly, all four lines bunch together with a slope of about −1. The reason is the jumps of density at layer boundaries: a piece that straddles a jump of height J errs by an amount of order J·h, whatever rule is applied to it, and there is always such a piece at every boundary. Nine boundaries, nine errors of order h, and the order of any rule falls to one.

The error without the answer: Runge’s rule

In the table above the error is easy to compute: the answer is known. In a real problem it is not, or there would be no point in computing. What helps is that the method’s error behaves predictably: if the order is p, then with half the step it decreases 2^p times. So two answers with steps h and h/2 differ mainly by the error, and from their difference it can be recovered.

This is Runge’s rule: compute twice and divide the difference by 2^p − 1, by 3 for trapezoids and by 15 for Simpson. The last column of the table above is exactly such an estimate of Simpson’s error, without knowing the exact answer: at 4 pieces it missed by almost a factor of two, and from 32 pieces on it matches the real error to within a per cent. The proviso in the claim matters: the rule works when h is already small enough for the error to behave as Ch^p. On a jump or on noisy data there is no such regime, and the estimate becomes a rough guess.

The estimate can not only be printed but also added to the answer: I_h/2 + (I_h/2 − I_h)/3 for trapezoids is exactly Simpson’s rule, and the error becomes of order h⁴. Repeating the trick with corrections of ever higher orders gives Romberg’s method, a table where each next column is more accurate than the previous. The trick itself, refining by two steps, was described in general form by Lewis Richardson in 1911; textbooks call it Richardson extrapolation.

An adaptive step: subdividing where it is hard

An equal step is wasteful. In the Earth the density inside a layer changes smoothly, and a few points would be enough there, while at each of the nine boundaries thousands are needed. An adaptive method applies Runge’s rule not to the whole interval but to each piece. It computes Simpson on a piece whole and on its two halves. If the answers agree with room to spare, the piece is done. If not, it splits it in half and repeats with each half, giving it half the tolerance.

Step through it on √x. Its derivative at zero is infinite, and the method goes ever deeper towards zero: the strip of points under the axis thickens on the left, while on the right a few pieces suffice. At a tolerance of 10⁻⁷ it needs 185 evaluations and errs by 8·10⁻¹⁰; Simpson with an equal step on the same 185 points errs by 3·10⁻⁵. Switch to the Earth: the method knows nothing of layers but finds all nine boundaries by itself, subdividing the interval near each one right to the bottom, like bisection looking for a root.

The number 15 in the stopping condition comes from Runge’s rule: Simpson has p = 4, and the error of the two halves is roughly 2⁴ − 1 = 15 times smaller than their discrepancy with the whole. The difference divided by 15 the method adds to the answer at once, like Richardson. The tolerance is halved on splitting so that the sum of the errors of all pieces stays within the original. On the Earth blindly it makes 1485 evaluations at a tolerance of 10⁻⁶ and gets an error of 1.3·10⁻⁹, almost a hundred thousand times less than Simpson with an equal step on 4097 points. The depth of subdivision is limited to fifty: on a jump the discrepancy never disappears, and without a limit the method would split a piece down to the size of rounding error.

The recursive adaptive Simpson was first published by G. F. Kuncir in 1962, as algorithm 103 in Communications of the ACM. Its present form, with halving and a correction of one fifteenth, was given to it by James Lyness in 1969.

A trap: a spike between the nodes

Switch the scene to the narrow spike: the parabola 1 + x² with a Gaussian spike of width 0.002 at the point 0.37. The method looks at the points 0, 1/4, 1/2, 3/4 and 1, sees a parabola on which Simpson’s rule is exact (the halves agree with the whole to the last digit) and stops after five evaluations with the answer 1.3333333. The correct answer is 1.3368782. Make the tolerance 10⁻¹² if you like: the answer is the same. Move the spike to 1/4 or 1/2 or make it wider, and the method will notice it and subdivide.

So no quadrature guarantees accuracy without assumptions about the function. A guarantee appears only if something is known about the function in advance: its derivative is bounded, there are no spikes narrower than some width, the places of jumps are known. Libraries ask exactly this: scipy.integrate.quad has a points parameter for the difficult places inside the interval, singularities and discontinuities. This is the very “every layer on its own” switch from the first scene.

Gauss nodes: twice the accuracy for the same evaluations

In all the rules above the nodes are equally spaced, and the only freedom is in the weights. A rule on n points with well-chosen weights is exact for polynomials of degree up to n − 1: the weights are chosen so that the integral matches on 1, x, …, xⁿ⁻¹. Gauss asked: what if the nodes themselves are moved too? Then there are twice as many free numbers, n nodes and n weights, and one can hope for exactness up to degree 2n − 1. The hope comes true if the nodes are put at the roots of a Legendre polynomial.

The program finds the roots of P_n by Newton’s method, starting from the cosines of an even grid of angles, and computes P_n itself by the three-term recurrence (k + 1)P_k+1 = (2k + 1)xP_k − kP_k−1. Three nodes, 0 and ±√(3/5) with weights 8/9 and 5/9, integrate x⁰, …, x⁵ exactly and first err on x⁶. PREM in each layer is a polynomial of degree at most three, and after multiplication by r² at most five. So three Gauss points per layer give the mass of the Earth exactly, and ten layers are thirty evaluations.

The example at the end of the program is taken from Gauss’s own 1815 memoir, “A new method of approximately computing integrals”. The integral of dx/ln x from 100,000 to 200,000 interested Gauss for a reason: it roughly equals the number of primes in this range. Gauss computed it with one, two, …, seven points by hand, and his seven-point answer, 8406.2431211, is correct to nine digits. In our program seven nodes give an error of 2.6·10⁻⁹, while Simpson on the same seven points errs by 7·10⁻³. There are actually 8392 primes between 100,000 and 200,000, but that is a question about primes, not about quadrature.

nodes   Gauss            error        Simpson on as many nodes
    3   8406.2367752    -6.346e-03      4.421e-01
    5   8406.2431171    -3.779e-06      3.399e-02
    7   8406.2431208    -2.550e-09      7.061e-03

The Gauss rule has an inconvenience: the nodes for n and for 2n points do not coincide, and when refining all the evaluations have to be thrown away. In 1964 Alexander Kronrod found how to add n + 1 more nodes to the n Gauss nodes so that the new rule is exact up to degree 3n + 1 and the old evaluations are put to use. The difference of the two answers is an error estimate almost for free. The QUADPACK library of 1983 rests on this pair, and through it scipy.integrate.quad: adaptive subdivision with a 21-point Kronrod rule on each piece.

When there is only a table of measurements

Everything above required the right to evaluate the function anywhere. With measurements this is impossible: the points stand where they stand, and neither an adaptive step nor Gauss nodes are available. Trapezoids and Simpson remain, and the question changes: not “which rule” but “are there enough measurements”.

An example is heating. How much heat a house loses over a winter is roughly proportional to degree-days: the integral of the difference between 18 °C and the outdoor temperature over all hours of the heating season. Below is Almaty from 15 October 2024 to 15 April 2025, by the hour. Over all 4368 hours it comes to 3036 degree-days. Now imagine an observer who records the temperature once a day and assumes it stayed that way until the next record.

A record at 14:00 undercounts almost 27% of the heat: by day it is warmer than the daily mean, and the observer sees the warmest hour every day. A record at 6 in the morning overestimates by 19%. At 9 in the morning or 7 in the evening the error is under a per cent, but that is luck: at those hours the temperature passes through its daily mean. Look at the line under the chart: trapezoids and Simpson on the same records give almost the same error. A high-order rule does not help, because the error here comes not from how the points are joined but from when they were taken: the step equals the period of the daily temperature cycle, and all records fall at the same phase of it.

Four records a day, at 5, 11, 17 and 23 hours (that is 0, 6, 12 and 18 UTC, the main observing times of weather stations), err by 0.7%, eight by 0.2%. This is not small by chance. The daily cycle repeats every 24 hours, and a sum over N equally spaced points over a period is exact for a constant and for any harmonic whose number is not divisible by N: the sum of e^2πikj/N over j from 0 to N − 1 is zero, as is the integral. Four records make no error on waves with periods of 24, 12 and 8 hours, and the daily cycle has almost all its amplitude in the first harmonics. For a smooth periodic function such a sum converges not as h² but faster than any power; Lloyd Trefethen and André Weideman devoted a 2014 survey to this, “The exponentially convergent trapezoidal rule”. Simpson on a periodic function is even worse than trapezoids: its weights of 4 and 2 spoil the equality of the points.

Runge’s rule works on a table too: compute the integral over all measurements and over every second one. Here that is 3036.22 and 3036.24 degree-days, so hourly measurements are more than enough. If the two answers differ by per cents, no rule will save you: more frequent measurements are needed. And if measurements come irregularly, with the sensor skipping hours, trapezoids work at any step: each trapezoid takes its own width. That is how numpy.trapezoid and scipy.integrate.simpson are built, which are given the measurement times along with the values.

Code: the Earth by layers

Here is the program behind the first scene. The density is the PREM polynomials, the exact answer comes from the antiderivatives of each layer, and then four rules blindly and by layers. The same three Gauss points give everything else seismologists know about the interior from one density profile.

moment of inertia I/(MR^2) = 0.3308
gravity: 9.82 m/s^2 at the surface, the most 10.69 at 3480 km from the centre
pressure in the centre: 364.1 GPa

The moment of inertia 0.3308MR² is less than the 0.4 of a uniform ball: the mass is concentrated towards the centre. This number is measured independently from the precession of the Earth’s axis, and the NASA fact sheet gives the same 0.3308. Gravity at depth is computed from the mass inside radius r, again an integral, and it is surprising: down to the core boundary at a depth of 2891 km it does not fall but grows to 10.69 m/s². The shell overhead pulls neither up nor down, and the mass underfoot decreases more slowly than the factor 1/r² grows: almost a third of the Earth’s mass is in the heavy core. The pressure at the centre is the weight of all the layers above it, the integral of ρ(r)g(r), where g is itself an integral: 364 gigapascals, three and a half million atmospheres. An integral inside an integral with the same three points per piece: inside a layer everything is smooth, and 64 pieces per layer are more than enough.

What it costs

The time of any quadrature is the number of function evaluations: summing costs almost nothing, while the function may be heavy, up to a whole simulation per point. The memory is O(1) for constant-step rules and a stack as deep as the subdivision for the adaptive one. So the main question is: how many evaluations are needed for accuracy ε.

On jumps and kinks all columns but the last lose their meaning: the order falls to one unless the interval is cut at the singular points. If you know where they are, cut; if not, the adaptive method will find them itself, at the price of extra evaluations.

In several dimensions the price grows fast. A grid of n points along each of d axes is n^d evaluations, and Simpson’s error in terms of the number of points N falls as N^−4/d. For the Monte Carlo method it falls as N^−1/2 in any dimension, and at d = 8 the two laws are equal. Beyond that random points beat any grid.

Where it works

• Scientific libraries. scipy.integrate.quad is a wrapper over QUADPACK: adaptive subdivision, the Gauss–Kronrod rule on pieces, default tolerances of about 1.5·10⁻⁸ and the points parameter for discontinuities.
• Tables of measurements. numpy.trapezoid and scipy.integrate.simpson take measurement times along with values: energy from power readings, distance from speed, volume from flow rate.
• The finite element method. The strength of a bridge or a wing is computed by dividing the structure into elements; the integrals over each element are taken not by formula but at a few Gauss points.
• Pharmacology. How much of a drug reached the blood is estimated by the area under the concentration curve (AUC) from samples taken at unequal intervals, by trapezoids.
• Construction. The heating degree-days of a region set how warm walls must be: the same sum as in the last scene, only over many years of data.

History

Trapezoids are older than the integral by almost two thousand years. On cuneiform clay tablets from 350–50 BCE astronomers computed how far Jupiter would move across the sky in 60 days as the area of a trapezoid under the graph of its daily speed, and found the moment when half the way was covered by dividing the trapezoid into two of equal area. The historian Mathieu Ossendrijver deciphered these calculations in 2016. Kepler proposed a parabola instead of a straight line in 1615 for the volume of wine barrels; Thomas Simpson published the general formula in 1743, more on which in the chapter on the integral of “Mathematics, the Queen of the Sciences”. Equal-step rules up to eleven points were derived by Roger Cotes; his papers were published in 1722, after the author’s death, and since then these rules have been called the Newton–Cotes rules.

Gauss’s memoir on free nodes appeared in 1815, and he derived them through continued fractions. The derivation through orthogonal polynomials, by which they are taught now, was given by Carl Jacobi in 1826. The rule of estimating the error by two steps was proposed by Carl Runge at the start of the twentieth century, Lewis Richardson turned it into a way of refining the answer in 1911, and Werner Romberg in 1955 into a table of repeated refinements. The sixties brought computer quadrature: Kuncir’s adaptive Simpson in 1962, Kronrod’s rules in 1964, and in 1983 Robert Piessens, Elise de Doncker-Kapenga, Christoph Überhuber and David Kahaner assembled the QUADPACK library, which still works inside SciPy.

And one story about how slowly knowledge spreads. In 1994 the journal Diabetes Care published a paper by the nutrition researcher Mary Tai: a way to compute the area under a blood sugar curve by dividing it into trapezoids and adding them up. It was the trapezoidal rule, and the journal soon printed letters with telling titles such as “Tai’s formula is the trapezoidal rule”. Tai replied that she had derived the method herself in 1981 and published it because colleagues were already using it and wanted something to cite. They still cite it: by 2025 the paper had more than five hundred citations.

Problems

In the first problem real Almaty hours with gaps: the integral has to be taken over unequal steps and carefully where the temperature crosses 18 °C. In the second, a planet of layers: the mass, moment of inertia and pressure at the centre, where one integral sits inside another.