Section4.2Numerical Integration
We want to build methods for approximating integrals. Two of the most common are the
Trapezoidal rule and Simpson's rule. Each of these integration approximation techniques
requires us to partition the domain into small bit and then to approximate the function
with either linear or curved polynomial functions. We will use polynomials since the
exact area can be worked out analytically.
We want to approximate \(\displaystyle \int_a^b f(x) dx\). One of the simplest ways is
to approximate the area under the function with a trapezoid. Recall from basic
geometry that area of a
trapezoid is \(A = \frac{1}{2} (b_1 + b_2) h\). In terms of the integration problem we
can do the following:
First parition \([a,b]\) into the set \(\{x_0=a, x_1, x_2, \ldots, x_{n-1},
x_n=b\}\).
On each part of the partition approximate the area with a trapezoid:
\begin{equation*}
A_j = \frac{1}{2} \left[ f(x_j) + f(x_{j-1}) \right]\left(
x_j - x_{j-1} \right)
\end{equation*}
Approximate the integral as
\begin{equation*}
\int_a^b f(x) dx \approx \sum_{j=1}^n A_j
\end{equation*}
Draw a picture depicting how the trapezoidal rule works.
The trapezoidal rule does a decent job approximating integrals, but ultimately you are using linear functions to approximate \(f(x)\) and the accuracy may suffer if the step size is too large or the function too non-linear. You likely notice that the trapezoidal rule will give an exact answer if you were to integrate a linear or constant function. An potentially better approach would be to get an integral that evaluates quadratic functions exactly. In order to do this we need to
evaluate the function at three points (not two like the trapezoidal rule). Let's to integrate a function \(f(x)\) on the interval \([a,b]\) by using the three points \( (a,f(a)), \, (m,f(m)),\) and \((b,f(b))\) where \(m =
\frac{a+b}{2}\) is the midpoint of the two boundary points. We want to find constants \(A_1, A_2, A_3\) such that the integral
\begin{equation*}
\int_a^b f(x) dx = A_1 f(a) + A_2 f\left( \frac{a+b}{2}\right) + A_3 f(b)
\end{equation*}
is exact for all constant, linear, and quadratic functions. This would guarantee that we have an exact method for all polynomials of order 2 or less but should serve as a decent approximation if the function is not quadratic.
To find the constants \(A_1, A_2, A_3\) we can write the following system of three equations
\begin{equation*}
\int_a^b 1 dx = b-a = A_1 + A_2 + A_3 \\
\int_a^b x dx = \frac{b^2-a^2}{2} = A_1 \cdot a + A_2 \cdot \left( \frac{a+b}{2} \right) + A_3 \cdot b \\
\int_a^b x^2 dx =\frac{b^3-a^3}{3} = A_1 \cdot a^2 + A_2 \cdot \left( \frac{a+b}{2} \right)^2 + A_3 \cdot b^2.
\end{equation*}
Solving the linear system gives
\begin{equation*}
A_1 = \frac{b-a}{6} \quad A_2 = \frac{4(b-a)}{6} \quad \text{and} \quad A_3 = \frac{b-a}{6}.
\end{equation*}
At this point we can see that the integral can be approximated as
\begin{equation*}
\int_a^b f(x) dx \approx \frac{b-a}{6} \left( f(a) + 4 f\left( \frac{a+b}{2} \right) + f(b) \right)
\end{equation*}
and the technique will give an exact answer for any polynomial of order 2 or below. To improve upon this idea we now examine the problem of partitioning the interval \([a,b]\) into small pieces and running this process on each piece.
Now we put the process explained above into a form that can be coded to approximate integral. We call this method Simpson's Rule after Thomas Simpson (1710-1761) who, by the way, was a basket weaver in his day job so he could pay the bills and keep doing math.
First parition \([a,b]\) into the set \(\{x_0=a, x_1, x_2, \ldots, x_{n-1},
x_n=b\}\).
On each part of the partitions approximate the area with a parabola:
\begin{equation*}
A_j = \frac{1}{6} \left[ f(x_j) + 4 f\left( \frac{x_j+x_{j-1}}{2} \right) +
f(x_{j-1}) \right]\left( x_j - x_{j-1} \right)
\end{equation*}
Approximate the integral as
\begin{equation*}
\int_a^b f(x) dx \approx \sum_{j=1}^n A_j
\end{equation*}
Draw a picture depicting how Simpson's rule works.
Write MATLAB functions that implement both the trapezoidal rule and Simpson's rule.
Your functions should accept, as input, an anonymous function handle, the bounds of
integration, and the number of interior points used in the integration.
Keep in mind that MATLAB deals with vectors and iteration in very nice ways. You
shouldn't need a for loop in your function.
Test both pieces on known integrals and approximate the order of the error based on
the mesh size.
Write MATLAB function that implements the trapezoidal rule,
but this time the only input is a list of \((x,y)\) ordered pairs representing the
integrand function. (Note: why aren't we doing Simpson's rule this way?)
Go to the http://apps.who.int/gho/data/?theme=home (The World Health Organization
Data Repository) and find a data set where it would make physical sense to integrate
(and that the integral would tell you something meaningful about the data). Use your
MATLAB code to take the integral and provide context and meaning associated with the
integral. You will need to download the data as a CSV or Excel file and import the
proper columns into Excel (Google “xlsread” to learn how MATLAB imports Excel
data). Your data set will obviously involve some noise, but one thought might be to
find a data set that shows a trend in time.
Consider problem 7 from Chapter 10 of the text. Instead of making a table, create a
loglog plot that shows the errors for the integral with different values of \(h\) (log
of \(h\) on the \(x\)-axis and log of the error on the \(y\)-axis).
Write a complete interpretation of the loglog plot. Then create a similar plot for
the integral
\begin{equation*}
\int_{-2}^2 e^{-x^2/2} dx.
\end{equation*}
(Coding Challenge) The four adjacent digits in the 1000-digit number that have
the greatest product are \(9 \times 9 \times 8 \times 9 = 5832\).
Write code to find the thirteen adjacent digits in the 1000-digit number that have the
greatest product. What is the value of this product?