Approximate an integral with a controlled error, integrals from easy to advanced lesson

The function $e^{-x^2}$ is smooth and positive on $[0, 1]$ , so the area $\int_0^1 e^{-x^2}\,dx$ is a definite number, but none of the usual functions gives an antiderivative to evaluate. Sampling the curve and adding weighted heights gives an estimate, and the real question is how many of its digits can be trusted. This page builds composite Simpson’s rule from parabolas through sample points. It then uses Simpson’s error bound, together with a bound on the fourth derivative, to prove that the integral equals $0.746824$ to six decimal places.

Key terms

No elementary antiderivative

Definition The elementary functions are the functions built from polynomials, roots, exponentials, logarithms, trigonometric functions and their inverses by finitely many arithmetic operations and compositions. A function has no elementary antiderivative when none of its antiderivatives is an elementary function. A theorem of Liouville shows that $e^{-x^2}$ is such a function, and so are $\frac{\sin x}{x}$ and $\frac{1}{\ln x}$ . Having no elementary antiderivative is not the same as having no antiderivative. If $f$ is continuous on an interval containing $0$ , the fundamental theorem of calculus makes $F(x) = \int_0^x f(t)\,dt$ an antiderivative of $f$ on that interval.

In plain words The area under $e^{-x^2}$ from $0$ to $1$ exists and is one definite number. What is missing is a closed formula for an antiderivative, so the number has to be computed another way. The substitution $u = -x^2$ does not help, because $du = -2x\,dx$ needs a factor $x$ that the integrand does not have.

Example With the factor $x$ present, the antiderivative is elementary,

\int x e^{-x^2}\,dx = -\tfrac12 e^{-x^2} + C.

Without it, the antiderivative is usually written with the error function $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2}\,dt$ , which gives

\int_0^1 e^{-x^2}\,dx = \frac{\sqrt{\pi}}{2} \operatorname{erf}(1).

The error function is defined by this very integral, so the formula names the number without computing it. Its digits still have to come from an approximation.

Numerical integration

Definition Numerical integration approximates a definite integral $\int_a^b f(x)\,dx$ by a weighted sum of function values

\sum_{j=0}^{n} w_j f(x_j)

at finitely many nodes $x_j$ in $[a, b]$ , with weights $w_j$ fixed by the rule. Riemann sums, the trapezoidal rule and Simpson’s rule are rules of this kind. A numerical value is an approximation until an error bound, or another proof, says how close it is.

In plain words Sample the integrand at a few points, weight the samples and add. Rules differ in where they sample and how they weight. A digit that stays the same as a rule is refined, or two rules that agree, is a useful hint but not a proof, because both can miss the same feature between the samples.

Example On $[0, 1]$ the function $e^{-x^2}$ decreases from $1$ to $e^{-1} \approx 0.3679$ , so the crudest rules, a single rectangle with the height at either end, already enclose the value,

e^{-1} < \int_0^1 e^{-x^2}\,dx < 1.

With two subintervals of width $\tfrac12$ , the trapezoidal rule joins the samples at $0$ , $\tfrac12$ and $1$ by straight lines and gives

\begin{aligned} T_2 &= \tfrac14\bigl(1 + 2e^{-1/4} + e^{-1}\bigr) \\[0.4em] &\approx 0.731370, \end{aligned}

which is too small by about $0.0155$ . The next card uses the same three samples with different weights and does far better.

Quadratic interpolant

Definition Given three points $(x_0, f_0)$ , $(x_1, f_1)$ and $(x_2, f_2)$ with $x_0 < x_1 < x_2$ , there is exactly one polynomial $p$ of degree at most $2$ whose graph passes through all three. It is called the quadratic interpolant of the points. When the nodes are equally spaced, with $x_1 - x_0 = x_2 - x_1 = h$ , its integral is

\begin{aligned} &\int_{x_0}^{x_2} p(x)\,dx \\[0.4em] &\quad = \frac{h}{3}\bigl(f_0 + 4f_1 + f_2\bigr). \end{aligned}

This is Simpson’s rule on a single pair of subintervals. To see the formula, put the middle node at $0$ with $t = x - x_1$ and write $p = A + Bt + Ct^2$ . The odd term $Bt$ integrates to $0$ over $[-h, h]$ , so the integral is $2Ah + \tfrac23 Ch^3$ . The interpolation conditions give $A = f_1$ and $2Ch^2 = f_0 - 2f_1 + f_2$ , and substituting these values gives the formula.

In plain words Over two neighboring subintervals, replace the curve by the parabola through its three sample points and take the exact area under the parabola. The weights $1$ , $4$ and $1$ come from integrating that parabola. They weight the samples in the area formula and do not change the heights themselves. The weights add to $6$ , so a constant height $c$ receives the correct area $2hc$ .

Example For $\int_0^1 e^{-x^2}\,dx$ with $h = \tfrac12$ , the samples at $0$ , $\tfrac12$ and $1$ are $1$ , $e^{-1/4}$ and $e^{-1}$ , and

\begin{aligned} S_2 &= \tfrac16\bigl(1 + 4e^{-1/4} + e^{-1}\bigr) \\[0.4em] &\approx 0.7471804. \end{aligned}

The integral is about $0.7468241$ , so this single parabola is off by about $3.6 \times 10^{-4}$ , more than $40$ times closer than the trapezoidal estimate from the same three samples.

Composite Simpson’s rule

Definition Let $f$ be defined on $[a, b]$ , let $n$ be an even positive integer, and set $h = (b - a)/n$ , $x_j = a + jh$ and $f_j = f(x_j)$ for $j = 0, 1, \dots, n$ . Composite Simpson’s rule applies the one-pair formula to $[x_0, x_2]$ , $[x_2, x_4]$ and so on up to $[x_{n-2}, x_n]$ , and adds the results,

\begin{aligned} S_n &= \frac{h}{3}\bigl(f_0 + 4f_1 + 2f_2 \\[0.4em] &\quad + 4f_3 + \dots + 2f_{n-2} \\[0.4em] &\quad + 4f_{n-1} + f_n\bigr). \end{aligned}

The endpoints get weight $1$ , the odd-numbered nodes get weight $4$ , and the interior even-numbered nodes get weight $2$ , because each of them is the shared end of two neighboring pairs. The pattern is written for $n \ge 4$ , and for $n = 2$ it is the one-pair formula. The number of subintervals must be even so that every subinterval belongs to exactly one pair. When $f$ is continuous on $[a, b]$ , the values $S_n$ approach $\int_a^b f(x)\,dx$ as $n$ grows through even numbers.

In plain words Cover the interval with a chain of parabola patches, one for each pair of neighboring subintervals, and add their areas. An odd number of subintervals would leave one subinterval without a partner.

Example Take the integral from the opening, whose integrand $e^{-x^2}$ is smooth on $[0, 1]$ , with $n = 4$ and $h = \tfrac14$ . The samples at $0$ , $\tfrac14$ , $\tfrac12$ , $\tfrac34$ and $1$ are about $1$ , $0.939413$ , $0.778801$ , $0.569783$ and $0.367879$ . The two odd-numbered samples add to about $1.509196$ , so

\begin{aligned} S_4 &\approx \tfrac{1}{12}\bigl(1 + 4(1.509196) \\[0.4em] &\quad + 2(0.778801) \\[0.4em] &\quad + 0.367879\bigr) \\[0.4em] &\approx 0.746855, \end{aligned}

an error of about $3.1 \times 10^{-5}$ . With $n = 100$ the same rule, carried out by a computer, gives $S_{100} \approx 0.7468241329$ , while a high-precision evaluation of the error function gives $\int_0^1 e^{-x^2}\,dx \approx 0.7468241328$ . The next card shows how to be sure of such digits without knowing the true value.

Simpson’s error bound

Definition Let $f$ have a continuous fourth derivative on $[a, b]$ , let $M$ be a number with $\lvert f^{(4)}(x) \rvert \le M$ for every $x$ in $[a, b]$ , and let $n$ be even with $h = (b - a)/n$ . Then composite Simpson’s rule satisfies

\begin{aligned} &\left\lvert \int_a^b f(x)\,dx - S_n \right\rvert \\[0.4em] &\quad \le \frac{(b - a)h^4 M}{180}. \end{aligned}

Since $h = (b - a)/n$ , the bound also equals $M(b - a)^5/(180 n^4)$ .

In plain words The bound is a guarantee that holds whenever its hypotheses do, and it needs no knowledge of the true value. It shrinks like $h^4$ , so halving $h$ divides it by $16$ . It describes the worst case, and the actual error is often much smaller. When the fourth derivative is not continuous on the whole interval, the theorem does not apply, and the error can be far larger than the formula would suggest.

Example Write $I = \int_0^1 e^{-x^2}\,dx$ . Here $b - a = 1$ , and the next card shows that $M = 76$ is a valid bound. With $n = 100$ ,

\begin{aligned} \lvert I - S_{100} \rvert &\le \frac{76}{180} \cdot 10^{-8} \\[0.4em] &< 4.23 \times 10^{-9}. \end{aligned}

So the integral lies between $S_{100} - 4.23 \times 10^{-9} \approx 0.7468241287$ and $S_{100} + 4.23 \times 10^{-9} \approx 0.7468241371$ . Every number in that interval rounds to $0.746824$ , so the integral equals $0.746824$ to six decimal places. The actual error is about $8 \times 10^{-11}$ , far below the guarantee.

Fourth derivative bound

Definition A fourth derivative bound for $f$ on $[a, b]$ is a number $M$ with $\lvert f^{(4)}(x) \rvert \le M$ for every $x$ in $[a, b]$ , not only at the sample points. Any valid $M$ may be used in Simpson’s error bound. When $f^{(4)}$ is continuous on $[a, b]$ , the smallest valid choice is the maximum of $\lvert f^{(4)} \rvert$ on $[a, b]$ , and a larger $M$ is still correct but calls for more subintervals.

In plain words Simpson’s rule integrates every cubic polynomial exactly. On a pair of subintervals centered at $0$ , the term $t^3$ has integral $0$ over $[-h, h]$ , and the rule also gives $\tfrac{h}{3}\bigl((-h)^3 + 0 + h^3\bigr) = 0$ . So the error depends only on how far $f$ is from being a cubic near each pair, and the fourth derivative, which is $0$ for every cubic, measures that.

Example For $f(x) = e^{-x^2}$ , differentiating four times gives

\begin{aligned} f^{(4)}(x) &= q(x)\,e^{-x^2}, \\[0.4em] q(x) &= 16x^4 - 48x^2 + 12. \end{aligned}

On $[0, 1]$ each power of $x$ is at most $1$ and $e^{-x^2} \le 1$ , so the triangle inequality gives $\lvert f^{(4)}(x) \rvert \le 16 + 48 + 12 = 76$ . This bound is quick to justify, and it is the one used above. A closer look gives the sharp bound $12$ . The fifth derivative $-8x(4x^4 - 20x^2 + 15)e^{-x^2}$ changes sign in $(0, 1)$ only at $x \approx 0.9586$ , so $f^{(4)}$ falls from $12$ at $x = 0$ to about $-7.42$ and then rises to $-20/e \approx -7.36$ at $x = 1$ . To guarantee an error below $10^{-6}$ , the bound $M = 76$ calls for $n = 26$ , while $M = 12$ needs only $n = 18$ .

Check yourself

Approximate $\int_0^2 x^4\,dx$ by Simpson’s rule with $n = 2$ , and compare the actual error with Simpson’s error bound.

Show the answer

Here $h = 1$ and the samples at $0$ , $1$ and $2$ are $0$ , $1$ and $16$ , so

\begin{aligned} S_2 &= \tfrac13(0 + 4 \cdot 1 + 16) \\[0.4em] &= \tfrac{20}{3}. \end{aligned}

The exact value is $\tfrac{32}{5}$ , so the error is $\tfrac{20}{3} - \tfrac{32}{5} = \tfrac{4}{15} \approx 0.267$ . The fourth derivative is the constant $24$ , so $M = 24$ and the bound is

\frac{2 \cdot 1^4 \cdot 24}{180} = \frac{4}{15}.

The bound is attained exactly here, so in general it cannot be lowered.

How many subintervals does composite Simpson’s rule need to guarantee an error below $10^{-6}$ for $\int_1^2 \frac{dx}{x}$ ?

Show the answer

For $f(x) = 1/x$ the fourth derivative is $f^{(4)}(x) = 24/x^5$ , which is continuous on $[1, 2]$ and largest in size at $x = 1$ , so $M = 24$ . With $b - a = 1$ and $h = 1/n$ the bound is

\frac{24}{180 n^4} = \frac{2}{15 n^4}.

It is below $10^{-6}$ when $n^4 > \tfrac{2}{15} \cdot 10^6 \approx 133333$ , that is, when $n > 19.1$ . The number $n$ must also be even, so $n = 20$ . The bound for $n = 20$ is about $8.3 \times 10^{-7}$ , while $n = 18$ only guarantees about $1.3 \times 10^{-6}$ . In fact $S_{20} \approx 0.6931474$ and $\ln 2 \approx 0.6931472$ , an actual error of about $1.9 \times 10^{-7}$ .

A classmate wants $\int_0^1 \sqrt{x}\,dx$ from Simpson’s rule, first with $n = 3$ and then with $n = 4$ together with the error bound. What goes wrong?

Show the answer

With $n = 3$ the three subintervals cannot be grouped into pairs, so composite Simpson’s rule is not defined. It needs an even number of subintervals. With $n = 4$ the rule applies. The samples at $0$ , $\tfrac14$ , $\tfrac12$ , $\tfrac34$ and $1$ are $0$ , $\tfrac12$ , $\tfrac{\sqrt2}{2}$ , $\tfrac{\sqrt3}{2}$ and $1$ , and with the weights $1, 4, 2, 4, 1$

\begin{aligned} S_4 &= \frac{3 + \sqrt{2} + 2\sqrt{3}}{12} \\[0.4em] &\approx 0.656526, \end{aligned}

compared with the exact value $\tfrac23 \approx 0.666667$ , an error of about $0.0101$ . The error bound cannot certify any of this. The fourth derivative $f^{(4)}(x) = -\tfrac{15}{16} x^{-7/2}$ is unbounded as $x \to 0^+$ , so it is not continuous on $[0, 1]$ and no finite $M$ exists. The error also shrinks slowly. Quadrupling $n$ to $16$ and then to $64$ gives errors of about $1.3 \times 10^{-3}$ and $1.6 \times 10^{-4}$ , a factor of about $8$ each time, while for an integrand that meets the hypotheses the bound shrinks by a factor of $4^4 = 256$ .

Approximate an integral with a controlled error

Key terms

No elementary antiderivative

Numerical integration

Quadratic interpolant

Composite Simpson’s rule

Simpson’s error bound

Fourth derivative bound

Check yourself

How this lesson connects

Used later in