Mathematics 134 – Intensive Calculus for Science 2 Lab Day 2 – “In search of a better numerical integration method” February 27, 2006

Background

Yesterday in class we discussed the LEF T (n), RIGHT (n), MID(n) and T RAP (n) methods for approximating definite integrals (the left-, right-, and midpoint Riemann sums were not new; the trapezoidal method was new). We have discussed the following patterns (some more than once!):

• If f is increasing on [a, b], then LEF T (n) gives an underestimate of ∫_{a}^{b} f(x) dx for all n. If f is decreasing on [a, b], then LEF T (n) gives an overerestimate of ∫_{a}^{b} f(x) dx for all n.

• If f is increasing on [a, b], then LEF T (n) gives an undererestimate of ∫_{a}^{b} f(x) dx for all n. If f is decreasing on [a, b], then LEF T (n) gives an overerestimate of ∫_{a}^{b} f(x) dx for all n.

• Whether T RAP (n) is an under- or over-estimate of ∫_{a}^{b} f(x) dx depends on the concavity of f: If f is concave up on [a, b], then T RAP (n) gives an overerestimate of ∫_{a}^{b} f(x) dx for all n. If f is concave down on [a, b], then T RAP (n) gives an undererestimate of ∫_{a}^{b} f(x) dx for all n.

• Whether MID(n) is an under- or over-estimate of ∫_{a}^{b} f(x) dx also depends on the concavity of f, but in the opposite way from T RAP (n): If f is concave up on [a, b], then MID(n) gives an undererestimate of ∫_{a}^{b} f(x) dx for all n. If f is concave down on [a, b], then MID(n) gives an overerestimate of ∫_{a}^{b} f(x) dx for all n. (Recall we saw this by replacing the midpoint rectangles by “tangent trapezoids” with the same areas.)

Today, we want to gather some data on these methods by looking at several examples. We will also see an even better method obtained by combining two of them in an appropriate way.

Maple Commands and Examples

The commands for finding the left, right, and midpoint sums are the same as in Lab 1. You will need to start by entering the command:

                          with(student);

to load the student package. The ones we will use in this time are just the ones that compute numerical approximations (not the ones that show the graphics). For instance,

             evalf(leftsum(x^2 - 3*x + 4, x=0..2, 5));
            evalf(rightsum(x^2 - 3*x + 4, x=0..2, 5));
            evalf(middlesum(x^2 - 3*x + 4, x=0..2, 5));

compute the LEF T (5), RIGHT (5), and MID(5) approximations to ∫_{0}^{2} x^{2} − 3x + 4 dx. There is a similar command for the trapezoidal rule. This does T RAP (5) for the same function as above:

                evalf(trapezoid(x^2 - 3*x + 4, x=0..2, 5));

As you can probably guess now, the format for all of these commands is: the command name, open parenthesis, the formula for the function f as a function of x, comma, x =, then the endpoints, separated by two periods, another comma, then the number n, followed by the close parenthesis, then the semicolon.

We will also need to be able to get exact values (or at least very close approximations) to our integrals. This is done in Maple by commands like this:

                        int(t^2 - 3*t + 4, t=0..2);
                    evalf(Int(t^2 - 3*t + 4, t=0..2));

Try these and look closely at the output. The first applies the FTC and gives the exact value. The second applies Maple’s “super-accurate” numerical methods to give a decimal approximation that is correct to 8 or 9 decimal places at least. (Note the capital I on the Int here – it’s important, but it’s slightly complicated to explain exactly what it means – i.e. “don’t ask” unless you really want to get a peak “under the hood” at what Maple actually does with your input commands(!).)

There will be some cases where Maple will not be able to find an antiderivative of the f you give it; in that case the output will be the same integral back again. For instance try

                          int(exp(x^3), x=0..1);

This means that Maple was unable to find an elementary antiderivative for the function f(x) = e^{x^{3}}, so it could not carry out the FTC to find the definite integral. (In fact this is an example where no elementary antiderivative exists.) To approximate that integral, we would use

                       evalf(Int(exp(x^3), x=0..1));

A Worked Example

Suppose we want to compare the T RAP (20) approximation to ∫_{0}^{1} e^{x^{3}} dx with the exact value. Here is a Maple session illustrating a good way to organize the computation. First assign the name f to the function (so we don’t need to retype the formula over and over):

                               f:=exp(x^3);

Then we compute the trapezoidal rule approximation, and assign it a name (trap20):

                  trap20:=evalf(trapezoid(f,x=0..1,20));

We take the results from Maple’s “super-accurate” numerical method as the exact value:

                exact:=evalf(Int(f,x=0..1));

Then we use Maple as a numerical calculator to compute the error for the approximation (exact value) - (approx. value):

                    err:=exact - trap20;

The error here should be −0.001696932. The size here indicates that the first three decimal places of the approximate value are correct (if we round).

Lab Problems

A) For each of the following integrals,

1) Compute an accurate numerical approximation using the evalf(Int(function, limits)); command as described above. We will treat this as our exact value – it’s the most accurate estimate we know!

2) Compute LEF T (n), RIGHT (n), MID(n), and T RAP (n) approximations for n = 20, 40, 80, 160, 320, and compute the error for each. For each integral, arrange your data into a table, giving the approximation and the error for the four methods and the five different n values. Format: Use one row in your table for each method, and put in columns for the values for each n and the associated error (with n increasing left to right).

Integrals:

1) ∫_{0}^{2} e^{−x^{2}/2} dx (enter the function as exp(-x^2/2))

2) ∫_{1}^{2} \frac{sin x}{x} dx

3) ∫_{0}^{1} \sqrt{1 + x^{4}} dx (enter the function as sqrt(1+x^4))

B) Now we want to look for some patterns in our data.

1) For each integral and each method separately, do you notice any consistent pattern when you compare the size of the error with a given n and with n twice as large (e.g. the error for MID(20) vs. the error for MID(40), or the error for T RAP (160) vs. the error for T RAP (320))? Is the pattern the same for all of the methods, or does it vary? NOTE: you may need to ignore the first entry for each method before the patterns become apparent!

2) Do you notice any consistent pattern when you compare the sizes of the errors for the four different methods on the same integral, with the same n? In particular, what is the approximate relation between the size of the errors for the TRAP and MID methods (for the same integral and the same n)?

C) One commonly-used better integration method is called Simpson’s Rule (Homer Simpson??) One way to write the formula for Simpson’s rule is:

SIMP (n) = \frac{2 · MID(n) + T RAP (n)}{3}

There is another command called simpson in the student package in Maple that uses this method to compute approximate values of integrals.

1) Try it on the examples from question A, and compare the sizes of the errors for Simpson’s Rule and the other methods for each n.

2) Why is Simpson’s Rule apparently more accurate? (Hint: Think about your answer to part 2 of question B).

Assignment

One lab write-up from each pair, due Wednesday, March 1.