Math 451: Computational Mathematics
3.6 Simpson's Rule
3.6.1 Newton-Cotes Formulas ¶
The Trapezoid Rule, or Newton-Cotes Formula of order 1, uses narrow trapezoids to approximate the area under the curve of a function like $f(t) = 3\cos(\pi t) + 5$.
However the Trapezoid Rule doesn't converge to the area quickly.
Instead let us use the midpoint of an interval to get three points and find the unique parabola that passes through those three points.
We approximate the area under the curve with an
interpolating polynomial
; in this case a quadratic polynomial. It is the unique quadratic polynomial that passes through
$$
A = (.5, f(.5)),\quad M = \left(\dfrac{.5 + 4}{2}, f\left(\dfrac{.5 + 4}{2}\right)\right),\ \text{and}\ B = (4, f(4)).
$$
This approximation is not very accurate. With two subintervals we obtain,
The visual representation of the approximation with four subintervals results in a much better approximation.
3.6.1 The Interpolating Polynomial ¶
There is a unique polynomial of degree $n-1$ that passes through $n$ distinct points on the plane. This is due to the fact that a polynomial of degree $n-1$ has $n$ coefficients that must be determined.
$$
y = c_{n-1}x^{n-1} + c_{n-2}x^{n-2} + \dots + c_1x + c_0
$$
To find a quadratic polynomial that approximates our function on an interval $[a,b]$, we will need three points on the curve of the function; typically we choose the midpoint $m = \dfrac{a+b}{2}$. The system of three equations and three unknowns can be written
$$
\begin{align*}
c_2a^2 + c_1a + c_0 &= f(a) \\
\\
c_2m^2 + c_1m + c_0 &= f(m) \\
\\
c_2b^2 + c_1b + c_0 &= f(c) \\
\\
\end{align*}
$$
In matrix form our system is written
$$
\begin{bmatrix} a^2 & a & 1 \\ m^2 & m & 1 \\ b^2 & b & 1 \end{bmatrix}\begin{bmatrix} c_2 \\ c_1 \\ c_0 \end{bmatrix} = \begin{bmatrix} f(a) \\ f(m) \\ f(b) \end{bmatrix}
$$
The matrix $A = \begin{bmatrix} a^2 & a & 1 \\ m^2 & m & 1 \\ b^2 & b & 1 \end{bmatrix}$ is called the
Van Der Monde Matrix
and the system has a unique solution as long as $a$, $b$ and $m$ are all have different values. We could solve this system for a complicated expression for $c_2$, $c_1$ and $c_0$ and then integrate the polynomial to obtain the area under the curve from $a$ to $b$. However there is an equally complicated formulation of the interpolating quadratic polynomial discovered by
Joseph-Louis Lagrange
,
$$
P(x) = f(a)\dfrac{(x-m)(x-b)}{(a-m)(a-b)} + f(m)\dfrac{(x-a)(x-b)}{(m-a)(m-b)} + f(c)\dfrac{(x-a)(x-m)}{(b-a)(b-m)}
$$
We recognize $P$ as a polynomial of degree two because each term is a polynomial of degree two. Moreover $P(a) - f(a) + 0 + 0$, while $P(m) = 0 + f(m) + 0$ and $P(b) = 0 + 0 + f(b)$. This is the correct formulation of the interpolating polynomial because we want to calculate the area between the curve and the horizontal axis. As we will see, the formula for the integrals simplifies to a simple formula for the area.
At this point we should play close attention to plus and minus signs. As $a \lt m \lt b$, we have that $m-b\lt 0$, $a-m\lt 0$, and $a-b\lt 0$. Switching two factors in the first term and one factor in the second term yields
$$
\begin{align*}
\text{Area} &= \dfrac{f(a)}{(m-a)(b-a)}\displaystyle\int_a^b (x-m)(x-b)\,dx \\
\\
&\qquad - \dfrac{f(m)}{(m-a)(b-m)}\displaystyle\int_a^b (x-a)(x-b)\,dx \\
\\
&\qquad + \dfrac{f(b)}{(b-a)(b-m)}\displaystyle\int_a^b (x-a)(x-m)\,dx
\end{align*}
$$
Now we compute first of these integrals using
integration by parts
twice.
$$
\begin{align*}
u &= x-m &\ dv &= x-b \\
du &= 1 &\ v &= \dfrac{(x-b)^2}{2} \\
du &= 0 &\ v &= \dfrac{(x-b)^3}{6} \\
\end{align*}
$$
The integral becomes
$$
\begin{align*}
\displaystyle\int_a^b (x-m)(x-b)\,dx &= \left[ (x-m)\dfrac{(x-b)^2}{2} - 1\dfrac{(x-b)^3}{6} \right]
_a^b \\
&= \left[ 0 - 0 - \left( \dfrac{(a-m)(a-b)^2}{2} - \dfrac{(a-b)^3}{6} \right) \right] \\
&= \dfrac{(m-a)(b-a)^2}{2} - \dfrac{(b-a)^3}{6} \\
\end{align*}
$$
If we recall that $m = \dfrac{a+b}{2}$, then we have that $m-a = \dfrac{a+b}{2} - a = \dfrac{b-a}{2}$. Substituting the value of our integral and our expression for $m$ into the first term results in
$$
\begin{align*}
\dfrac{f(a)}{(m-a)(b-a)}\displaystyle\int_a^b (x-m)(x-b)\,dx &= \dfrac{2f(a)}{(b-a)^2}\left[ \dfrac{(b-a)^3}{4} - \dfrac{(b-a)^3}{6} \right] \\
&= f(a)(b-a)\left[ \dfrac{1}{2} - \dfrac{1}{3} \right] \\
&= f(a)\dfrac{b-a}{6} \\
\end{align*}
$$
The second integral evaluates to
$$
\begin{align*}
\displaystyle\int_a^b (x-a)(x-b)\,dx &= \left[\dfrac{(x-a)(x-b)^2}{2} - \dfrac{(x-b)^3}{6} \right]_a^b \\
&= \left[ 0 - 0 - \left( 0 - \dfrac{(a-b)^3}{6} \right) \right] \\
\end{align*}
$$
As with the first term $b-m = b - \dfrac{a+b}{2} = \dfrac{b-a}{2}$. Hence the second term becomes
$$
\begin{align*}
-\dfrac{f(m)}{(m-a)(b-m)}\displaystyle\int_a^b (x-a)(x-b)\,dx &= -\dfrac{4f(m)}{(b-a)^2}\left[ -\dfrac{(b-a)^3}{6} \right] \\
&= 4f(m)\dfrac{b-a}{6} \\
\end{align*}
$$
The third integral in
$$
\begin{align*}
\displaystyle\int_a^b (x-a)(x-m)\,dx &= \left[\dfrac{(x-m)(x-a)^2}{2} - \dfrac{(x-a)^3}{6} \right]_a^b \\
&= \left[ \dfrac{(b-m)(b-a)^2}{2} - \frac{(b-a)^3}{6} - \left( 0 - 0 \right) \right] \\
&= \dfrac{(b-a)^3}{4} - \frac{(b-a)^3}{6} \\
&= \dfrac{(b-a)^3}{12}
\end{align*}
$$
Hence the third term becomes
$$
\begin{align*}
\dfrac{f(b)}{(b-a)(b-m)}\displaystyle\int_a^b (x-a)(x-m)\,dx &= \dfrac{2f(b)}{(b-a)^2}\cdot\dfrac{(b-a)^3}{12} \\
&= f(b)\dfrac{b-a}{6}
\end{align*}
$$
Therefore our total area is given by
$$
\begin{align*}
\text{Area} &= \displaystyle\int_a^b P(x)\,dx \\
&= f(a)\dfrac{b-a}{6} + 4f(m)\dfrac{b-a}{6} + f(b)\dfrac{b-a}{6} \\
\\
&= \dfrac{b-a}{6}\left[ f(a) + 4f(m) + f(b) \right] \\
\\
&= \dfrac{b-a}{6}\left[ f(a) + 4f\left(\dfrac{a+b}{2}\right) + f(b)\,\right] \\
\end{align*}
$$
3.6.3 Simpson's Rule of Order 2 ¶
Simpson's rule consists of approximating the area under a curve of a function on an interval $[a,b]$ by calculating the area of the interpolating quadratic polynomial on the interval. One obtains a more accurate approximation by dividing the interval $[a,b]$ into smaller subintervals and using Simpson's Rule to compute the areas for each the interpolating polynomials for each subinterval.
If we partition the interval $[a,b]$ into $n$ subintervals of equal length, then each of the subintervals of length $h = \dfrac{b-a}{n}$. This gives us $n+1$ points
$$
\begin{align*}
a = x_0 &< x_1 = x_0+h \\
&< x_2 = x_0+2h \\
&< x_3 = x_0 + 3h \\
&< \dots \\
&< x_{n-1} = x_0 + (n-1)h \\
&< x_n = x_0 + nh = a + n\dfrac{b-a}{n} = b
\end{align*}
$$
The approximate area using Simpson's Rule for each subinterval is given by
$$
\begin{align*}
\text{Area}_1 &= \dfrac{x_1-x_0}{6}\left[ f(x_0) + 4f\left(\dfrac{x_0 + x_1}{2}\right) + f(x_1) \,\right] \\
\\
\text{Area}_2 &= \dfrac{h}{6}\left[ f(x_1) + 4f\left(\dfrac{x_1 + x_2}{2}\right) + f(x_2) \,\right] \\
\\
\ddots \\
\\
\text{Area}_k &= \dfrac{h}{6}\left[ f(x_{k-1}) + 4f\left(\dfrac{x_{k-1} + x_k}{2}\right) + f(x_k)\,\right] \\
\\
\ddots \\
\\
\text{Area}_n &= \dfrac{h}{6}\left[ f(x_{n-1}) + 4f\left(\dfrac{x_{n-1} + x_n}{2}\right) + f(x_n)\,\right]
\end{align*}
$$
Adding these $n$ areas to get the total area yields
$$
\begin{align*}
\text{Area} &= \dfrac{h}{6} \left[ f(x_0) + 4f\left(\dfrac{x_0 + x_1}{2}\right) + f(x_1) \right. \\
&\qquad + f(x_1) + 4f\left(\dfrac{x_1 + x_2}{2}\right) + f(x_2) \\
&\qquad\ddots \\
&\qquad + f(x_{k-1}) + 4f\left(\dfrac{x_{k-1} + x_k}{2}\right) + f(x_k) \\
&\qquad\ddots \\
&\qquad + \left. f(x_{n-1}) + 4f\left(\dfrac{x_{n-1} + x_n}{2}\right) + f(x_n)\,\right]
\end{align*}
$$
Notice that the right endpoint of a subinterval is the left endpoint of the next subinterval. This allows us to write the equation for the area
$$
\begin{align*}
\text{Area} &= \dfrac{h}{6} \left[ f(x_0) + 4f\left(\dfrac{x_0+x_1}{2}\right) + 2f(x_1) + 4f\left(\dfrac{x_1+x_2}{2}\right) + 2f(x_2) + \dots \right. \\
\\
&\qquad\left. \dots +\ 2f(x_{n-1}) + 4f\left(\dfrac{x_{n-1}+x_n}{2}\right) + f(x_n)\,\right]
\end{align*}
$$
This formula can be vectorized to compute the area for $n$ subintervals in one step.
Exercise 3.6.1 ¶
1. Write function simpson.m with the following declaration ¶
function [area, loopcount] = simpson(f, a, b, tolerance, maximum_loop_count)
simpson computes the area under the curve of an integrable function
f, on an interval [a,b] using the general Simpson's Rule to the specified tolerance.
Inputs:
f a handle to an integrable function
a left end point of interval of integration
b right end point of interval of integration
tolerance accuracy requested; use a default tolerance of sqrt(eps)
maximum_loop_count the maximum number of loops before exiting the algorithm
use a default maximum_loop_count of 1024
Outputs:
area the area between the curve and the x-axis
n the number of used to compute the Simpson's rule approximation
You should continue computing the general Simpson's rule with an increasing number of subintervals until either the maximum_loop_count is reached or until the difference between successive estimates for the area is less than the tolerance.
Start with 4 (four) subintervals and double the number of subintervals until your algorithm terminates.
2. Write function quadrature.m with the following declaration ¶
function area = quadrature(f, a, b, n, tolerance)
quadrature computes the area under the curve of an integrable function
f, on an interval [a,b] using the adaptive Simpson's rule for a maximum
recursion depth n, or to the specified tolerance.
Inputs:
f a handle to an integrable function
a left end point of interval of integration
b right end point of interval of integration
n the maximum depth of recursion
tolerance the required acceptable error in the computation of the area
Outputs:
area the area between the curve and the x-axis
Remember that you will only need Simpson's rule for
one
interval $[a,b]$ in the recursion step.
$$
\dfrac{b-a}{6}\left[ f(a) + 4f\left(\dfrac{a+b}{2}\right) + f(b)\,\right]
$$
3. Compare the results of simpson.m and quadrature.m for the following functions ¶
- $.5 + \cos(x)$ on the interval $\left[0,\ \frac{\pi}{4}\right]$
- $.5 + \sin(x)$ on the interval $\left[-\pi,\ 0\right]$
- $\tan\left(x-\frac{\pi}{4}\right)$ on the interval $\left[0,\ \frac{\pi}{2}\right]$
- $\text{acos}(x) - .5$ on the interval $[-1,\ 1]$
- $\text{asin}(x - .5)$ on the interval $[-.5,\ 1.5]$
- $\text{atan}(x - .5)$ on the interval $[-64,\ 64]$
- $@(x)\ \text{polyval}\left(\begin{bmatrix} -12 & 4 & 15 & -5 & -3 & 1 \end{bmatrix}, x\right)$ on the interval $[-1.5 1.5]$
- $\log(\text{abs}(x))$ on the interval $[0.5,\ 32]$
- $\text{abs}(x) - \pi$ on the interval $[-8,\ 8]$
You should compare these functions in a MATLAB script or function that you write so that I can just run and see your results.
Your use of this self-initiated mediated course material is subject to our Creative Commons License .
Creative Commons Attribution-NonCommercial-ShareAlike 4.0
Attribution
You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
Noncommercial
You may not use the material for commercial purposes.
Share Alike
You are free to share, copy and redistribute the material in any medium or format. If you adapt, remix, transform, or build upon the material, you must distribute your contributions under the
same license
as the original.