Skip to main content

Numerical Analysis A look at rootfinding, quadrature and approximation using analytic and computing tools

Peter Staab

Section 5.3 Composite Newton-Cotes Quadrature

In the previous section, we developed a class of numerical integration techniques called Newton-Cotes Quadrature formulas. In this case we break down an integral into smaller subintervals and then apply the formulas from the previous section.
🔗
If we let \(x_{i} = a+i h\) where \(h=(b-a)/n\text{,}\) for \(n\) some positive integer, we can write
\begin{equation*} \begin{aligned}\int_{a}^{b} f(x) \, dx\amp = \int_{x_0}^{x_1}f(x) \, dx +\int_{x_1}^{x_2}f(x) \, dx + \cdots + \int_{x_{n-1}}^{x_n}f(x) \, dx \\\amp \approx \widehat{I}_{1}(f) + \widehat{I}_{2}(f) + \cdots + \widehat{I}_{n}(f)\end{aligned} \end{equation*}
where \(\widehat{I}_{k}\) is some quadrature formula on \([x_{k-1},x_{k}]\text{.}\) Note, in this section, the subscript does not determine the quadrature formula (like trapezoid or Simpson’s). Instead it only determines the interval on which the formula is applied.
🔗

Subsection 5.3.1 Composite Trapezoid

Recall that the Trapezoid Rule with its error term from Theorem 5.2.9 is
\begin{equation*} \begin{aligned}\int_{a}^{b} f(x) \, dx\amp = \frac{b-a}{2}\bigl(f(a) + f(b) \bigr) - \frac{(b-a)^{3}}{12}f''(\xi)\end{aligned} \end{equation*}
and applying this to the interval \([x_{j},x_{j+1}]\) is
\begin{equation*} \begin{aligned}\int_{x_j}^{x_{j+1}}f(x) \, dx\amp = \frac{x_{j+1}-x_{j}}{2}\bigl(f(x_{j}) + f(x_{j+1}) \bigr) - \frac{(x_{j+1}-x_{j})^{3}}{12}f''(\xi_{j})\\ \amp \qquad \text{and assuming that the points are equally-space, then $x_{j+1}-x_{j}=h$ and $\xi_{j} \in [x_{j},x_{j+1}]$}\\ \amp = \frac{h}{2}\bigl(f(x_{j}) + f(x_{j+1}) \bigr) - \frac{h^{3}}{12}f''(\xi_{j})\end{aligned} \end{equation*}
🔗
So applying this on \(n\) subintervals of length \(h\) is
\begin{align} \int_{a}^{b} (x) \,\amp = \int_{x_0}^{x_1}f(x) \, dx +\int_{x_1}^{x_2}f(x) \, dx + \cdots + \int_{x_{n-1}}^{x_n}f(x) \, dx, \notag\\ \amp \underbrace{\frac{h}{2} \bigl( (f(x_0) + f(x_1)) + (f(x_1) + f(x_2)) + (f(x_2)+ f(x_3)) + (f(x_3)+f(x_4)) + \cdots }_{\text{Trapezoid Rule}}\notag\\ \amp \qquad \underbrace{+ (f(x_{n-2}) + f(x_{n-1})) + (f(x_{n-1}) + f(x_n)) \bigr)}_{\text{Trapezoid Rule}}- \underbrace{\frac{h^{3}}{12} \sum_{j=1}^n f''(\xi_j) }_{\text{Error}}\tag{5.3.1} \end{align}
🔗
And this can be written more compactly.
🔗
The Composite Trapezoid Rule on \([a,b]\) for \(n\) equally-spaced subintervals is designated \(T_{n}(f)\) and is
\begin{equation} T_{n}(f) = \frac{h}{2}\biggl( f(x_{0}) + 2 \sum_{j=1}^{n-1}f(x_{j}) + f(x_{n}) \biggr),\tag{5.3.2} \end{equation}
where \(h=(b-a)/n\) and \(x_{j} = a+jh\text{.}\)
🔗
Let’s look at the error formula a bit more carefully. Let \(c_{1}\) and \(c_{2}\) be defined to be:
\begin{equation*} \begin{aligned}f''(c_{1})\amp = \max_{x \in [a,b]}f''(x),\amp f''(c_{2})\amp = \min_{x \in [a,b]}f''(x) ,\end{aligned} \end{equation*}
therefore
\begin{equation*} \begin{aligned}f''(c_{2}) \leq f''(\xi_{i}) \leq f''(c_{1}) ,\end{aligned} \end{equation*}
for each \(i\text{.}\) If we sum over \(i\) from 1 to \(n\text{,}\) we get:
\begin{equation*} \begin{aligned}n f''(c_{2}) \leq \sum_{i=1}^{n} f''(\xi_{i}) \leq n f''(c_{1}),\end{aligned} \end{equation*}
🔗
By the Intermediate Value Theorem there exists a number \(\xi \in [a,b]\) such that
\begin{equation*} \begin{aligned}n f''(\xi) = \sum_{i=1}^{n} f''(\xi_{i}).\end{aligned} \end{equation*}
🔗
The error at the end of (5.3.2) can now be written:
\begin{equation*} -\frac{h^{3}}{12}\sum_{j=1}^{n} f(\xi_{i}) = -\frac{h^{3}}{12}\cdot n f''(\xi) = -\frac{(b-a)h^{2}}{12}f''(\xi), \end{equation*}
Where \(h=(b-a)/n\) has been used. Therefore,
\begin{equation*} \begin{aligned}\int_{a}^{b} f(x) \, dx\amp = \frac{h}{2}\biggl( f(x_{0}) + 2 \sum_{i=1}^{n-1}f(x_{i}) + f(x_{n}) \biggr) - \frac{(b-a)h^{2}}{12}f''(\xi)\end{aligned} \end{equation*}
🔗
The error of the Composite trapezoid rule is therefore:
\begin{equation} \int_{a}^{b} f(x) \, dx - T_{n}(f) = - \frac{(b-a)h^{2}}{12}f''(\xi),\tag{5.3.3} \end{equation}
for some \(\xi \in [a,b]\text{.}\)
🔗

Example 5.3.1.

Find \(\int_{0}^{1} e^{-x^2}\,dx\) using the composite trapezoid rule with \(n=10\text{.}\) Find a bound on the error.
🔗
Solution.
In this case \(h=(b-a)/n=1/10\text{.}\)
\begin{equation*} \begin{aligned}\int_{0}^{1} e^{-x^2}\,dx\amp = \frac{1}{20}\biggl( f(0) + 2 \sum_{i=1}^{9} f(i/10) + f(1) \biggr) \\\amp \approx 0.7462107961\end{aligned} \end{equation*}
🔗
To find a bound on the error, we use (5.3.3) and need a bound on the second derivative of \(f(x)\text{.}\) It can be shown that \(f''(x) = (4x^{2}-2)e^{-x^2}\) and this reaches a maximum (in absolute value) at \(x=0\text{,}\) so \(|f''(x)| \leq 2\text{.}\) Using (5.3.3),
\begin{equation*} \begin{aligned}\biggl|- \frac{(b-a)h^{2}}{12}f''(\xi) \biggr| \leq \frac{1(1/10)}{12}\max_{x \in [0,1]}f''(\xi) \leq \frac{1}{120}(2) = \frac{1}{60}.\end{aligned} \end{equation*}
🔗
Recall that this is a theoretical error bound and the actual error is about \(6 \times 10^{-4}\text{.}\)
🔗
🔗
🔗
🔗

Subsection 5.3.2 Rate of Convergence

(5.3.3) shows that the error is \(O(h^{2})\text{.}\) The following example will also show this is true numerically.
🔗

Example 5.3.2.

Show numerically that the error of composite trapezoid rule is \(O(h^{2})\) by the following. Find \(T_{n}(f)\) for \(n=1,2,4,8,\ldots,256\) where \(f(x) = \cos x\) on \([0,\pi/2]\text{.}\) Build a table of errors and show that \(e_{2h}/e_{h}\) appears to converge to 4.
🔗
Solution.
First, we will find \(T_{n}(f)\) for \(n=1,2,4,8,\ldots,256\text{.}\) To demonstrate the first few, for \(n=1\text{,}\) \(h/2 = (b-a)/2=\pi/4\text{,}\) so
\begin{equation*} \begin{aligned}T_{1}(f)\amp = \frac{\pi}{4}(f(0) + f(\pi/2)) \approx 0.78539\end{aligned} \end{equation*}
and when \(n=2\text{,}\) \(h/2=(b-a)/4=\pi/8\text{,}\) so
\begin{equation*} \begin{aligned}T_{2}(f)\amp = \frac{\pi}{8}\bigl( f(0) + 2 f(\pi/4) + f(\pi)\bigr) \approx 0.948059\end{aligned} \end{equation*}
and so on. The error is given by \(e_{h}=I(f)-T_{n}(f)=1-T_{n}(f)\) and finally a ratio \(e_{2h}/e_{h}\) is calculated. The following table lists all these values for \(n=1,2,4, \ldots, 256\text{.}\)
🔗
\(n\) \(h\) \(T_{n}(f)\) \(e_{h}=I(f)-T_{n}(f)\) \(e_{2h}/e_{h}\)
\(1\) \(\pi/2\) \(0.7853981635\) \(0.2146018365\)
\(2\) \(\pi/4\) \(0.9480594490\) \(0.0519405510\) \(4.1316819400\)
\(4\) \(\pi/8\) \(0.9871158012\) \(0.0128841988\) \(4.0313372840\)
\(8\) \(\pi/16\) \(0.9967851722\) \(0.0032148278\) \(4.0077415030\)
\(16\) \(\pi/32\) \(0.9991966798\) \(0.0008033202\) \(4.0019257580\)
\(32\) \(\pi/64\) \(0.9997991945\) \(0.0002008055\) \(4.0004890300\)
\(64\) \(\pi/128\) \(0.9999498004\) \(0.0000501996\) \(4.0001414350\)
\(128\) \(\pi/256\) \(0.9999874508\) \(0.0000125492\) \(4.0002231220\)
\(256\) \(\pi/512\) \(0.9999968633\) \(0.0000031367\) \(4.0007651350\)
An \(O(h^{2})\) method will have \(e_{2h}/e_{h} \rightarrow 4\) because since
\begin{equation*} \begin{aligned}e_{h}\amp = c h^{2}\amp e_{2h}\amp = c (2h)^{2}\\ \amp \text{finding the ratio}\frac{e_{2h}}{e_{h}}\\ \amp = \frac{4ch^{2}}{ch^{2}}= 4 .\end{aligned} \end{equation*}
🔗
🔗
🔗
We will show later in this chapter that knowing that the trapezoid rule is \(O(h^{2})\) will allow us to create an accelerated quadrature method called Romberg integration.
🔗
🔗

Subsection 5.3.3 Composite Simpson’s Rule

Since Simpson’s Rule has 2 subintervals, when we extend it to a composite Simpson’s rule, we need an even number of subintervals for \(n\text{,}\) so let \(n=2m\text{,}\) and
\begin{equation*} \begin{aligned}h\amp = \frac{b-a}{n}= \frac{b-a}{2m}\\ x_{i}\amp = a + ih\end{aligned} \end{equation*}
for \(i = 0, 1, 2, \ldots, n\text{.}\)
🔗
If \(\widehat{I}_{k}(f)\) represents Simpson’s Rule on \([x_{2k},x_{2k+1}]\text{,}\) then
\begin{equation*} \begin{aligned}\int_{a}^{b} f(x) \,\amp = \int_{x_0}^{x_2}f(x) \, dx +\int_{x_2}^{x_4}f(x) \, dx + \cdots + \int_{x_{n-2}}^{x_n}f(x) \, dx \\\amp \approx \widehat{I}_{0} (f) + \widehat{I}_{1} (f) + \cdots + \widehat{I}_{m} (f) \\\amp = \underbrace{\frac{h}{3} \bigl( f(x_0) + 4 f(x_1) + f(x_2) + f(x_2) + 4f(x_3) + f(x_4) +}_{\text{Simpson's rule}}\\\amp \qquad \underbrace{\cdots+ f(x_{n-1}) + 4 f(x_{n-1}) + f(x_n) \bigr)}_{\text{Simpson's rule}}\\\amp \qquad - \underbrace{ \sum_{j=1}^n \frac{(x_{2j}-x_{2j-2})^{5}}{2880}f''''(\xi_j) }_{\text{Error}}\end{aligned} \end{equation*}
Composite Simpson’s Rule on the interval \([a,b]\) is given by
\begin{equation} S_{n}(f) = \frac{h}{3}\biggl( f(x_{0}) + 4 \sum_{j=1}^{m} f(x_{2j-1}) + 2 \sum_{j=1}^{m-1}f(x_{2j}) + f(x_{2m}) \biggr)\tag{5.3.4} \end{equation}
for \(n\) even and \(h=(b-a)/n\) and \(x_{i}=a+ih\text{.}\)The error associated with Simpson’s rule is
\begin{equation} \int_{a}^{b} f(x) \,dx - S_{n} (f) = - \frac{(b-a)h^{4}}{180}f^{(4)}(\xi) \tag{5.3.5} \end{equation}
for some number \(\xi \in [a,b]\text{.}\)
🔗
The error analysis above for Simpson’s rule is similar to that found in the composite trapezoid rule above.
🔗

Example 5.3.3.

Find \(\int_{0}^{1} e^{-x^2}\, dx\) using composite Simpson’s rule and \(m=5\text{.}\) (That is, \(n=10\)). Also, find the theoretical error bound using \(|f^{(4)}(x)| \leq 10\) on \([0,1]\text{.}\)
🔗
First, note that the nodes are \(x_{j}=jh\) for \(j=0,1,2, \ldots, 10\) and \(h=1/10\text{.}\) Let \(f(x) = e^{-x^2}\text{.}\)
\begin{equation*} \begin{aligned}I_{k}(f)\amp = \frac{1}{30}\bigl(f(0) + 4 f(0.1) + 2 f(0.2) + 4 f(0.3) + 2 f(0.4) + 4 f(0.5) \\\amp \qquad \qquad + 2 f(0.6) + 4 f(0.7) + 2 f(0.8) + 4 f(0.9) + f(1) \bigr) \\\amp \approx 0.7468249483.\end{aligned} \end{equation*}
🔗
The error can be found by using (5.3.5),
\begin{equation*} \begin{aligned}\biggl\vert \int_{a}^{b} e^{-x^2}\,dx - S_{10}(f) \biggr\vert\amp = \biggl\vert - \frac{(1-0) (1/10)^{4}}{180}f^{(4)}(\xi) \biggr\vert \\ \leq 5.555 \times 10^{-6}\end{aligned} \end{equation*}
so the approximate value of the integral given above is guaranteed accurate to 5 decimal places.
🔗
🔗
The methods given here in this section can be extended to the other Newton-Cotes rules like the 3/8-rule and Boole’s rule. However, they are generally not used. It is often simpler to stay with the composite trapezoid or Simpson’s rule and use more subintervals to decrease the error. Also as we will see at the end of this chapter an accelerated method using composite trapezoid method will often generate more accurate results with fewer function evaluations.
🔗
🔗
🔗
PrevTopNext