Subsections

• Purpose
• Background
• Maple commands
• Error term
• Exercises

Numerical Integration

Purpose

The purpose of this lab is to give you some experience with using the trapezoidal rule and Simpson's rule to approximate integrals.

Background

The trapezoidal rule and Simpson's rule are used for approximating area under a curve or the definite integral

$$\int_{a}^{b} f(x) \, dx$$

Both methods start by dividing the interval $[a,b]$ into $n$ subintervals of equal length by choosing a partition

$$a = x_0 < x_1 < x_2 < \ldots < x_n = b$$

satisfying

$$x_i = a + i h, \mbox{ for $i=0, \ldots, n$}$$

where

$$h = \frac{b-a}{n}$$

is the length of each subinterval. For the trapezoidal rule, the integral over each subinterval is approximated by the area of a trapezoid. This gives the following approximation to the integral

$$T_n = \frac{h}{2} \left( f(x_0) + 2 f(x_1) + 2 f(x_2) + \ldots + 2 f(x_{n-1}) + f(x_n) \right)$$

For Simpson's rule, the function is approximated by a parabola over pairs of subintervals. When the areas under the parabolas are computed and summed up, the result is the following approximation.

$$S_n = \frac{h}{3} \left( f(x_0) + 4 f(x_1) + 2 f(x_2) + \ldots + 4 f(x_{n-1}) + f(x_n) \right)$$

Maple commands

The commands for the trapezoidal rule and Simpson's rule are in the student package.

>with(student);

The following example will use the function

>f:=x->x^2*exp(x);

This computes the integral of the function from 0 to 2.

>int(f(x),x=0..2);

Using the evalf command provides a decimal approximation.

>evalf(int(f(x),x=0..2));

The command for using the trapezoidal rule is trapezoid. The syntax is very similar to that of the int command. the last argument specifies the number of subintervals to use. In the command below, the number of subintervals is set to 10, but you should experiment with increasing or decreasing this number. Note that Maple writes out the sum and doesn't evaluate it to a number.

>trapezoid(f(x),x=0..2,10);

Putting an evalf command on the outside computes the trapezoidal approximation.

>evalf(trapezoid(f(x),x=0..2,10));

The command for Simpson's rule is very similar.

>simpson(f(x),x=0..2,10);
>evalf(simpson(f(x),x=0..2,10));

Error term

There is an error term associated with the trapezoidal rule that can be u sed to estimate the error. More precisely, we have

$$\int_{a}^{b} f(x) \, dx = T_n + E_n$$

where

$$E_n = - \frac{(b-a)^3 f''(c)}{12n^2}$$

for some value $c$ between $a$ and $b$ that maximizes the second derivative in absolute value. Solving the error formula for $n$ guarantees a number of subintervals such that the error term $\vert E_n\vert$ is less than some desired tolerance $\varepsilon$. This gives:

$$n > \sqrt{\frac{(b-a)^3 M}{12 \varepsilon}}$$

The way to think about this result is that it gives a value for $n$ which guarantees that the error of the trapezoidal rule is less than the tolerance $\varepsilon$. It is generally a very conservative result.

Similarly, the number of subintervals for the simpson rule approximation to guarantee an error smaller than $\varepsilon$ is

$$n > \sqrt[4]{\frac{(b-a)^5 M}{180 \varepsilon}}$$

Exercises

1. For the function $f(x)=x \sin(x)+4$, use the trapezoidal rule formula to approximate the area under $f(x)$ over the interval $[0,4]$ ising $n=8$. Verify your answer using the trapezoid command in Maple. Repeat this exercise with simpson's rule and Maple's simpson command.

2. For the function $$f(x)=\frac{x}{1+x^3}$$ over the interval $[0,5]$, complete the following steps.

   (i)

   By using Maple's int and possibly evalf commands, find a good approximation to the integral of the function over the given interval.

   (ii)

   Find the minimum number of subintervals $n$ necessary to approximate the value of the definite integral with error no greater than $0.001$. Then do the same with Simpson's rule. Which is more accurate and why?

   (iii)

   Use the error estimate $E_n$ for the trapezoidal rule to find a value for $n$, the number of subintervals, that ensures $T_n$ is within $0.001$ of the answer in part (i). How close is this answer to the number of subintervals you found in part (ii)?
   (Hint: The part of the error term that is most difficult is $f''(c)$. To find the maximum value of the second derivative of $f$ over the interval $[a,b]$, an approximation based on the plot should suffice. To plot the second derivative of $f$ over the given interval, type the following command in Maple:

   plot(abs(diff(f(x),x,x)),x=0..5);
   

   Then, using your mouse, click on the part of the plot that appears to be the maximum. When you click on a plot, $x, y$ coordinates will appear in the upper left hand corner of your Maple window and the $y$ coordinate is used for the max.)

   (a)

   (iv)] Repeat part (iii) using Simpson's rule. (Hint: To maximize the fourth derivative of $f$ over the given interval, type the following command in Maple:

   plot(abs(diff(f(x),x,x,x,x)),x=0..5);