Section 5.2 The Basics of Numerical Integration or Quadrature
We now turn to a more-important aspect of numerical calculus, numerical integration or quadrature. It is not difficult to find a function that doesn’t have an antiderivative that can written in terms of our standard set of functions. For example, we will use the integral
\begin{equation}
\int_{0}^{a} e^{-x^2}\, dx\tag{5.2.1}
\end{equation}
The basic idea behind quadrature is that for a continuous function \(f\) on an interval \([a,b]\text{,}\) we wish to find
\begin{equation*}
I(f) = \int_{a}^{b} f(x) \, dx.
\end{equation*}
Most quadrature formulas take on the form:
\begin{equation}
I(f) \approx I_{n} (f) = \sum_{i=0}^{n} w_{i} f(x_{i})\tag{5.2.2}
\end{equation}
where \(x_{i}\) are called the quadrature points and \(w_{i}\) are called the weights. The various methods we will see in this section and subsequent sections will differ based on finding different values of the nodes and weights.
This may look similar to that of a Riemann sum from Calculus. You may recall that the Riemann sum (for a finite number of rectangles) is an approximation to the area under the curve or the definite integral. Typically, if \(w_{i}=h=x_{i+1}-x_{i}\text{,}\) that is the points used in evaluation are equally spaced, then (5.2.2) is the Riemann sum.
Subsection 5.2.1 Newton-Cotes Formulas
The first class of formulas for quadrature that we will see are called the Newton-Cotes Formulas. In this case, we will use interpolating polynomials (usually using equally-spaced points) to approximate the function, then integrate the polynomial.
First, write the interpolating polynomial in terms of the Lagrange polynomials, \(L_{n,i}(x)\text{,}\)
\begin{equation*}
P_{n} (x) = \sum_{i=0}^{n} f(x_{i}) L_{n,i}(x).
\end{equation*}
Let \(I_{n}(f)\) be the definite integral of \(P_{n}(x)\) on \([a,b]\) or
\begin{equation*}
\begin{aligned}I_{n}(f)\amp = \int_{a}^{b} \sum_{i=0}^{n} f(x_{i}) L_{n,i}(x) \, dx, \\\amp = \sum_{i=0}^{n} f(x_{i}) \int_{a}^{b} L_{n,i}(x) \, dx.
\amp \text{This can be written as} \\ I_{n}(f)\amp = \sum_{i=0}^{n} w_{i} f(x_{i}),\end{aligned}
\end{equation*}
where
\begin{equation*}
w_{i} = \int_{a}^{b} L_{n,i}(x) \, dx
\end{equation*}
We now look at specific cases for \(n\text{,}\) the order of the interpolating polynomial.
Subsection 5.2.2 \(n=1\) case
In this case, there are two points, \(x_{0}=a\) and \(x_{1}=b\text{.}\)
\begin{equation*}
\begin{aligned}L_{1,0}(x)\amp = \frac{b-x}{b-a},\amp L_{1,1}(x)\amp = \frac{x-a}{b-a},\end{aligned}
\end{equation*}
therefore the weights are:
\begin{equation*}
\begin{aligned}w_{0}\amp = \int_{a}^{b} \frac{b-x}{b-a}\, dx\amp w_{1}\amp = \int_{a}^{b} \frac{x-a}{b-a}\, dx.
\\\end{aligned}
\end{equation*}
Although these integrals are not difficult to do, they are easier if we change variables. Let \(x = a+ h t\) where \(h=b-a\text{,}\) then these integrals become:
\begin{equation*}
\begin{aligned}w_{0}\amp = h\int_{0}^{1} (1-t) \, dt\amp w_{1}\amp = h\int_{0}^{1} t \, dt \\\amp = \frac{h}{2}\amp w_{1}\amp = \frac{h}{2}\end{aligned}
\end{equation*}
Thus the quadrature formula is
\begin{equation*}
\begin{aligned}I(f) \approx I_{1}(f)\amp = \frac{h}{2}\bigl( f(a) + f(b) \bigr).\end{aligned}
\end{equation*}
This formula is called the trapezoid rule, because the integral is equivalent to the area of the trapezoid formed by the points \((0,f(a))\text{,}\) \((a,f(a))\text{,}\) \((b,f(b))\text{,}\) and \((0,f(b))\) as seen in the plot below.
The trapezoid is the shaded area in the plot above and approximates the integral of \(f\) from \(a\) to \(b\text{.}\)
The Trapezoid Rule is an approximation to \(\int_{a}^{b}f(x) \, dx\text{,}\) given by
\begin{equation}
I_{1}(f) = \frac{h}{2}\bigl(f(a)+f(b)\bigr).\tag{5.2.3}
\end{equation}
Subsection 5.2.3 \(n=2\) case
In this case, we will use a quadratic function (parabola) that interpolates the data and then will be integrated. To do this, we use the following 3 points to determine the interpolating polynomial:
\begin{equation*}
\begin{aligned}x_{0}\amp = a,\amp x_{1}\amp = \frac{a+b}{2},\amp x_{2}\amp = b,\end{aligned}
\end{equation*}
and we will define \(h = (b-a)/2\) to be half the total distance between \(a\) and \(b\text{.}\) The Lagrange polynomials through these points are:
\begin{equation*}
\begin{aligned}L_{2,0}(x)\amp = \frac{(x-x_{1})(x-x_{2})}{2h^{2}},\\ L_{2,1}(x)\amp = \frac{(x-x_{0})(x-x_{2})}{-h^{2}}, \\ L_{2,2}(x)\amp = \frac{(x-x_{0})(x-x_{1})}{2h^{2}},\end{aligned}
\end{equation*}
and the weights are:
\begin{equation*}
\begin{aligned}w_{0}\amp = \int_{a}^{b} L_{2,0}(x) \, dx,\amp w_{1}\amp = \int_{a}^{b} L_{2,1}(x) \, dx,\amp w_{2}\amp = \int_{a}^{b} L_{2,2}(x) \, dx.\end{aligned}
\end{equation*}
If we define \(x=a+th\text{,}\) these can be written:
\begin{equation*}
\begin{aligned} w_{0}\amp = \frac{h}{2}\int_{0}^{2} (t-1)(t-2) \, dt = \frac{h}{3},\\ w_{1}\amp = -h \int_{0}^{2} t(t-2) \, dt = \frac{4h}{3}, \\ w_{2}\amp = \frac{h}{2}\int_{0}^{2} t(t-1) \, dt = \frac{h}{3}, \\\end{aligned}
\end{equation*}
so the quadrature formula in this case using the above points and weights and substituting into (5.2.2) can be written as follows.
The quadrature formula that uses three equally-spaced points is given by
\begin{equation*}
I_{2}(f) = \frac{h}{3}\biggl( f(a) + 4 f\biggl(\frac{a+b}{2}\biggr) + f(b) \biggr)
\end{equation*}
and is called Simpson’s Rule.
We now see a example from the start of this section that uses the trapezoid and Simpson’s rule to evaluate an integral.
For comparison, we will later show that
\begin{equation*}
\begin{aligned}\int_{0}^{1} e^{-x^2}\, dx \approx 0.7468241328\end{aligned}
\end{equation*}
which shows that Simpson’s Rule is correct to 3 decimal places and the trapezoid rule is correct to only 1 digit.
Subsection 5.2.4 Other Newton-Cotes
Other Newton-Cotes formula can be found by taking more points. In the \(n=3\) case, there are 4 points and 4 weights and it can be shown that
\begin{equation*}
\begin{aligned}w_{0}\amp = \frac{3h}{8},\amp w_{1}\amp = \frac{9h}{8}, \\ w_{1}\amp =\frac{9h}{8},\amp w_{3}\amp = \frac{3h}{8},\end{aligned}
\end{equation*}
\begin{equation*}
I_{3}(f) = \frac{3h}{8}\bigl( f(x_{0}) + 3 f(x_{1}) + 3 f(x_{2}) + f(x_{3}) ,\bigr)
\end{equation*}
which is called the 3/8-rule.
Lastly, when \(n=4\text{,}\) one can show that the weights are
\begin{equation*}
\begin{aligned}w_{0} =w_{4}\amp = \frac{14h}{45},\amp w_{1}= w_{3}\amp = \frac{64h}{45},\amp w_{2}\amp = \frac{24 h}{45},\end{aligned}
\end{equation*}
and using \(h=(b-a)/4\) with \(x_{i}=a+ih\) in (eq:quad:basic) results in Boole’s Rule for quadrature
\begin{equation*}
I_{4}(f) = \frac{h}{45}\bigl( 14 f(x_{0}) + 64 f(x_{1}) + 24 f(x_{2}) + 64 f(x_{3}) + 14 f(x_{4}) \bigr).
\end{equation*}
Example 5.2.3.
In this case, \(h=1/3\) and
\begin{equation*}
I(f) = \frac{1}{8}\bigl( e^{0} +3 e^{-1/9}+3 e^{-4/9}+ e^{-1}\bigr) \approx 0.7469923196
\end{equation*}
which is accurate to nearly 4 decimal places and with Boole’s rule, \(h=1/4\) and
\begin{equation*}
I_{4}(f) = \frac{1}{4\cdot 45}\bigl(14 e^{0} +64 e^{-1/16}+24 e^{-1/4}+ 64e^{-9/16}+14e^{-1}\bigr) \approx 0.7468337099
\end{equation*}
which is accurate to 5 decimal places.
Subsection 5.2.5 Error Formulas for Newton-Cotes
In this section, we will find error formulas for the Newton-Cotes rules derived above. That is, we will be able to determine theoretically how close the quadrature formula (approximation) is to the actual integral.
Definition 5.2.4.
The degree of precision or accuracy of a quadrature rule \(I_{n}(f)\) is the positive integer \(m\) such that
\begin{equation*}
\begin{aligned}I(p)\amp = I_{n}(p) \qquad \text{for every polynomial $p$ of degree $\leq m$}\\ I(p)\amp \neq I_{n}(p) \qquad \text{for some polynomial $p$ of degree $m+1$}\\\end{aligned}
\end{equation*}
In short the degree of precision of a quadrature formula is the highest degree polynomial that is integrated exactly. Due to the linearly of integration (and hence quadrature), it is sufficient to integrate only powers of \(x\) to determine the degree of precision.
Example 5.2.5.
If a rule integrates 1, \(x\) and \(x^{2}\) exactly, but not \(x^{3}\text{,}\) then the degree of precision is 2.
Example 5.2.6.
Find the degree of precision for Simpson’s Rule.
We build a table of the integrals of \(x^{n}\) and Simpson’s Rule,
\begin{equation*}
I_{2}(f) = \frac{h}{3}\bigl( f(a) + 4 f((a+b)/2) + f(b) \bigr)
\end{equation*}
| \(n\) | \(f(x)\) | \(\int_{a}^{b} f(x) \, dx\) | \(I_{2}(f)\) |
| 0 | 1 | \(b-a\) | \(\frac{b-a}{2 \cdot 3}(1 + 4 + 1) =b-a\) |
| 1 | \(x\) | \(\frac{b^{2}-a^{2}}{2}\) | \(\frac{b-a}{2 \cdot 3}( a + 2(a+b) + b)= \frac{b^{2}-a^{2}}{2}\) |
| 2 | \(x^{2}\) | \(\frac{b^{3}-a^{3}}{3}\) | \(\frac{b-a}{6}(a^{2} + 4(a+b)^{2}/4 + b^{2})=\frac{b^{3}-a^{3}}{3}\) |
| 3 | \(x^{3}\) | \(\frac{b^{4}-a^{4}}{4}\) | \(\frac{b-a}{6}\bigl(a^{3} + 4 \frac{(a+b)^{3}}{8}+ b^{3}\bigr)= \frac{b^{4}-a^{4}}{4}\) |
| 4 | \(x^{4}\) | \(\frac{b^{5}-a^{5}}{5}\) | \(\frac{b-a}{6}\bigl(a^{4} + 4 \frac{(a+b)^{4}}{16}+ b^{4}\bigr) \neq \frac{b^{5}-a^{5}}{5}\) |
Since Simpson’s rule integrates \(x^{n}\) exactly for \(n=0,1,2,3\) but not \(x^{4}\text{,}\) the degree of precision is 3.
The next example shows that we can use the degree of precision to find the weights. This is an alternative to using the interpolation polynomial to find the \(w_{i}\text{.}\)
Example 5.2.7.
Find \(w_{0}\) and \(w_{1}\) such that the quadrature rule on \([a,b]\) for \(\int_{a}^{b} f(x) \, dx\) given by
\begin{equation*}
\widetilde{I}_{1}(f) = w_{0} f(a+h) + w_{1} f(a+2h)
\end{equation*}
has as high a degree of precision as possible and \(h=(b-a)/3\text{.}\)
Solution.
First, recall that
\begin{equation*}
\int_{a}^{b} x^{n} \, dx = \frac{b^{n+1}-a^{n+1}}{n+1}\qquad \text{for $n=0,1,2,3,\ldots$}
\end{equation*}
Therefore we need to apply this quadrature formula for \(f(x) = x^{n}\text{,}\) \(n\) as large as possible.
-
If \(n=0\text{,}\) \(f(x)=1\text{,}\) so we apply this to the given quadrature rule.\begin{equation*} b-a = w_{0} + w_{1} \end{equation*}
-
If \(n=1\text{,}\) then \(f(x)=x\) and applying this to the quadrature rule,\begin{equation*} \frac{b^{2}-a^{2}}{2} = w_{0} (a+h) + w_{1} (a+ 2h) \end{equation*}Solving the two equations above and using the definition of \(h\text{,}\) we get:\begin{equation*} w_{0} = w_{1} = \frac{b-a}{2} \end{equation*}
So the quadrature rule becomes:
\begin{equation*}
\widetilde{I}_{1}(f) = \frac{b-a}{2}\bigl( f(a+h) + f(a+2h) \bigr)
\end{equation*}
If we apply this formula to \(f(x)=x^{2}\text{,}\) we get
\begin{equation*}
\widetilde{I}_{1}(x^{2}) = \frac{b-a}{18}(5a^{2}+8ab+5b^{2})
\end{equation*}
which is not \((b^{3}-a^{3})/3\text{,}\) so this quadrature rule has degree of precision 1.
This is an example of a Newton-Cotes open quadrature formula because the points do not include the endpoints.
An open quadrature formula does not include the endpoints. This can be useful for calculating improper integrals. For example, \(\int_{0}^{1} (1/\sqrt{x}) \, dx\) exists even though the integrand evaluated at \(x=0\) does not.
Next, we examine the error associated with quadrature formulas, but first need the following theorem.
Theorem 5.2.8. Weighted Mean Value Theorem.
If \(f\) is continuous on \([a,b]\text{,}\) \(g\) is integrable on \([a,b]\) and \(g(x)\) does not change sign on \([a,b]\text{,}\) then there exists a number \(\xi \in [a,b]\) such that
\begin{equation*}
\int_{a}^{b} f(x) g(x) \, dx = f(\xi) \int_{a}^{b} g(x) \, dx .
\end{equation*}
Proof.
Suppose \(g(x) \geq 0\) on \([a,b]\text{.}\) Since \(f\) is continuous on \([a,b]\text{,}\) it reaches a min and max on the interval, or \(m \leq f(x) \leq M\) so
\begin{equation*}
m g(x) \leq f(x) g(x) \leq M g(x),
\end{equation*}
for all \(x \in [a,b]\text{.}\) Thus,
\begin{equation*}
\begin{aligned}m \int_{a}^{b} g(x) \, dx \leq \int_{a}^{b} f(x) g(x) \, dx \leq M \int_{a}^{b} g(x) \, dx .\end{aligned}
\end{equation*}
If \(\int_{a}^{b} g(x) \,dx = 0\) then \(\int_{a}^{b} f(x) g(x) \,dx\) also equals zero since we assumed that \(g(x) \geq 0\text{,}\) so any \(\xi \in [a,b]\) satisfies the theorem. Otherwise, \(\int_{a}^{b} g(x) \, dx \gt 0\text{.}\) Therefore
\begin{equation*}
m \leq \frac{\int_{a}^{b} f(x) g(x) \, dx}{\int_{a}^{b} g(x) \, dx}\leq M.
\end{equation*}
Applying the intermediate value theorem, there exists a number \(\xi\text{,}\) such that
\begin{equation*}
f(\xi) = \frac{\int_{a}^{b} f(x) g(x) \, dx}{\int_{a}^{b} g(x) \, dx}.
\end{equation*}
If \(g(x)<0\text{,}\) then a similar argument can be made, thus the theorem holds true.
Theorem 5.2.9. Error Formula for Trapezoid Rule.
If \(f \in C^{2}[a,b]\text{,}\) and \(I_{1}(f)\) is the trapezoid rule given in (5.2.3), then there exists an \(\hat{\xi}\in[a,b]\) such that
\begin{equation*}
\int_{a}^{b} f(x) \, dx = I_{1}(f) - \frac{(b-a)^{3}}{12}f''(\hat{\xi}) .
\end{equation*}
Proof.
\begin{equation*}
f(x) = P_{1} (x) + \frac{f''(\xi(x))}{2}(x-a)(x-b) .
\end{equation*}
If we integrate then
\begin{equation*}
\int_{a}^{b} f(x) \, dx = I_{1}(f) + \int_{a}^{b} \frac{f''(\xi(x))}{2}(x-a)(x-b) \, dx,
\end{equation*}
where \(\int_{a}^{b} P_{1}(x) \, dx = I_{1}(f)\) as shown above. Since \((x-a)(x-b)\) does not change sign, we can apply the weighted mean value theorem applied to the integral in the last term, there is an \(\hat{\xi}\) such that
\begin{equation*}
\begin{aligned}\int_{a}^{b} f(x) \, dx\amp = I_{1}(f) + \frac{f''(\hat{\xi})}{2}\int_{a}^{b} (x-a)(x-b) \, dx, \\\amp = I_{1}(f) - \frac{f''(\hat{\xi})}{12}(b-a)^{3} .\end{aligned}
\end{equation*}
This theorem is quite useful for finding error bounds when using the trapezoid rule to approximate a definite integral. It also gives the degree of precision for the trapezoid rule. If we integral \(f(x)=x^{n}\text{,}\) then the error given by the second term in Theorem Theorem 5.2.9 which is proportional to \(f''(\xi)\) is 0 for \(x^{0}\) and \(x^{1}\text{,}\) but not \(x^{2}\text{.}\) Therefore the degree of precision is 1.
Theorem 5.2.10. Error Formula for Simpson’s Rule.
If \(f \in C^{4}[a,b]\) and \(I_{2}(f)\) is the formula for Simpson’s rule applied to \(f\) on \([a,b]\text{,}\) then there exists a \(\hat{\xi}\in [a,b]\) such that
\begin{equation*}
\int_{a}^{b} f(x) \, dx = I_{2}(f) - \frac{(b-a)^{5}}{2880}f^{(4)}(\hat{\xi})
\end{equation*}
Using an argument similar to that that above about the trapezoid rule, this shows that Simpson’s Rule has degree of precision of 3. The proof of this is more difficult than that of the error associated with the trapezoid rule, so is skipped. See Isaacson:1966fj for more information.
Example 5.2.11.
Use the Error Formula for Simpson’s Rule to find a theoretical error bound of \(\int_{0}^{1} e^{-x^2}\, dx\text{.}\)
First, one can show that if \(f(x) = e^{-x^2}\) then
\begin{equation*}
f''''(x) = 12 e^{-x^2}-48x^{2} e^{-x^2}+ 16 x^{4} e^{-x^2}
\end{equation*}
The maximum of this function on the interval \([0,1]\) occurs when \(x=0\) and is 12, so
\begin{equation*}
|f''''(\xi)| \leq 12.
\end{equation*}
Thus the error, \(E_{2}(x)\) defined as
\begin{equation*}
\begin{aligned}E_{2}(x)\amp = \biggl|\int_{0}^{1} e^{-x^2}\, dx - I_{2}(f)\biggr|, \\\amp \leq \frac{(b-a)^{5}}{2880}|f''''(\xi)|, \\\amp = \frac{12}{2880}= 0.004285.\end{aligned}
\end{equation*}
This means that the value for \(\int_{0}^{1} e^{-x^2}\, dx\) in Example 5.2.2 is accurate to 0.004285.
The next theorem shows that the error formula for any closed Newton-Cotes quadrature rules takes on a certain form and thus the degree of precision can be derived.
Theorem 5.2.12. Error Formula for Closed Newton-Cotes Quadrature.
If \(I_{n}(f)\) denotes the Closed Newton-Cotes Quadrature Rule with \(n+1\) nodes, then
-
if \(n\) is even and \(f \in C^{n+2}[a,b]\text{,}\) then there exist constants \(c \in \R\) and \(\xi \in [a,b]\) such that\begin{equation*} I(f) = I_{n}(f) - c(b-a)^{n+3}f^{(n+2)}(\xi) \end{equation*}and the degree of precision of \(I_{n}(f)\) is \(n+1\text{.}\)
-
if \(n\) is odd and \(f \in C^{n+1}[a,b]\text{,}\) then there exists constants \(c' \in \R\) and \(\xi' \in [a,b]\) such that\begin{equation*} I(f) = I_{n}(f) - c'(b-a)^{n+2}f^{(n+1)}(\xi') \end{equation*}and the degree of precision is \(n\text{.}\)
The proof of this theorem is not given but can be found in Bradie:2006qy.
Example 5.2.13.
Use the Error formula for Closed Newton-Cotes quadrature rules above to determine the degree of precision of Boole’s rule.
In the case of Boole’s rule, \(n=4\text{,}\) thus the exists constants \(c\) and \(\xi\) such that
\begin{equation*}
I(f) = I_{4}(f) -c(b-a)^{7} f^{(6)}(\xi),
\end{equation*}
where \(I_{4}(f)\) is Boole’s rule. The error is proportional to \(f^{(6)}(\xi)\text{,}\) thus the degree of precision is 5.
Subsection 5.2.6 Exercises and Projects
Project 5.2.1.
In Example 5.2.7, we derived an open Newton-Cotes formula for quadrature. An alternative way to derive these is to do a similar analysis of quadrature that we did above for the closed Newton-Cotes formulas.
Investigate the Open Newton-Coates formula and provide examples of applying the formula for various functions. Note that these can be helpful for analyzing type II improper integrals (see a standard calculus book).
Project 5.2.2.
Consider an integral of the form:
\begin{equation*}
\int_{-\infty}^b f(x) \, dx
\end{equation*}
or
\begin{equation*}
\int_{a}^{\infty} f(x) \, dx
\end{equation*}
The way to handle these integrals is to do a transformation (substitution) that maps from \((-\infty,b]\) or \([a,\infty)\) to a finite interval, say \([c,d]\) and then perform a standard quadrature technique on the transformed integral.
Investigate some possible transformations, apply some standard techniques and find error formulas for the new quadrature formula.
