> restart;

>

Numerical Integration

Note: While reading this worksheet make sure that you re-execute each command by placing the cursor on the command line and pressing Enter. You should do this with every worksheet, but it is even more important with this one. If you don't re-execute commands defining the integration rules LEFT, RIGHT, TRAP, MID, SIMP, Maple will not recognize them and you won't be able to solve your homework problems. You have to re-execute the commands in the same Maple session in which you do your homework.

You have learned how to estimate definite integrals using left and right Riemann sums as well as three other methods: the trapezoid, midpoint and Simpson's rules. In certain situations it is possible to tell if your estimates are too high or too low without actually knowing the exact value of your integral. For example, if the function is concave up on the interval of integration, the trapezoid rule overestimates and the midpoint rule underestimates. In this worksheet we will explore some other aspects of these approximations. Refer to Sec.7.5, 7.6 of you textbook for background.

>

Numerical Integration Rules Using Maple

Maple has all the above mentioned rules of numerical integration programmed in. They are available in the " student " package. However, the rules in the " student " package are defined slightly differently than in our textbook. To be consistent with the text, we shall program Maple to calculate approximations corresponding to each of the rules the way our book does. Although you are not expected to learn Maple programming at this time, you may want to look at the structure of our little programs below. In these programs we define procedures in Maple for calculating the left and right Riemann sums, trapezoid approximation, midpoint approximation and Simpson's approximation, for any given function f, interval [a,b], and the number of subdivisions n. " proc " stands for procedure. In parentheses are the variables f,a,b,n which we will have to specify each time we apply any of the procedures. The name " local " preceeding h and i tells Maple that these variables are to be used as indicated in the procedure, and not to be confused with variables of the same name which may appear elsewhere.

> LEFT:= proc (f,a,b,n) local h,i; h:=(b-a)/n; evalf(sum(f(a+i*h)*h,i=0..(n-1))); end:

> RIGHT:= proc (f,a,b,n) local h,i; h:=(b-a)/n; evalf(sum(f(a+i*h)*h,i=1..n)); end:

> TRAP:= proc (f,a,b,n) evalf((LEFT(f,a,b,n)+RIGHT(f,a,b,n))/2); end:

> MID:= proc (f,a,b,n) local h,i; h:=(b-a)/n; evalf(sum(f(a+(i+1/2)*h)*h, i=0..(n-1))); end:

> SIMP:= proc (f,a,b,n) evalf((2*MID(f,a,b,n)+TRAP(f,a,b,n))/3); end:

Notice how these procedures reflect the mathematical definitions of the corresponding numerical integration rules. In order to apply each of the procedures defined above, you need to give your function a name , say:

> w:=x->x^2;

To calculate the corresponding numerical approximation of an integral, for example , for, say, n =50 divisions, you use the following syntax:

> LEFT(w,0,1,50);

> SIMP(w,0,1,50);

>

and so on. You have to give your function a name ; the syntax LEFT(x^2,0,1,50) will not work. (Do you see why by looking at the programs?)

Once we have the rules of numerical integration programmed in, we can begin conducting numerical experimentation with Maple.

The Errors in Midpoint and Trapezoid Approximations

Let's choose an interesting function whose antiderivative cannot be found in simple terms. For example, . Let's define the function in Maple:

> f:=x->sin(x^2);

Suppose we want to find the integral . For the purpose of error analysis, we ask Maple to find the exact value of the integral:

> A:=Int(f(x),x=0..1); ExactA:=value(A);

The exact value is in terms of some function "Fresnel" that we don't know. Let's use the command "evalf", which by default will give us an approximation of the exact value to ten decimal places. Let it be our reference "exact" value.

> exactA:=evalf(ExactA);

>

Let's explore the error of approximation for our integral A using the trapezoid and the midpoint rules for n=25 and n=50 divisions. Let's label the results of applying the trapezoid and midpoint rules by T1, T2, M1, M2, respectively.

> T1:=TRAP(f,0,1,25);

> T2:=TRAP(f,0,1,50);

> M1:=MID(f,0,1,25);

> M2:=MID(f,0,1,50);

Recall that the error of an approximation is defined as Error = Exact value- Approximate value . Let's calculate the errors corresponding to the approximations that we obtained above. We label the errors as ErrorT1 , ErrorT2 , etc.

> ErrorT1:=exactA-T1; ErrorT2:=exactA-T2;

> ErrorM1:=exactA-M1; ErrorM2:=exactA-M2;

>

Observe that the midpoint errors are of the opposite sign than the corresponding trapezoid errors, and the magnitude of each midpoint error is, roughly, half of the magnitude of the corresponding trapezoid error. Indeed, let's calculate the corresponding ratios. More precisely, let's calculate the ratios of absolute values of errors as we are only interested in ratios of magnitudes of errors. The syntax for the absolute value is " abs(ErrorM1) " , etc.

> evalf(abs(ErrorM1)/abs(ErrorT1)); evalf(abs(ErrorM2)/abs(ErrorT2));

>

More numerical experimentation, with integrals of different functions over different intervals, confirms the above observation. This leads to the following rule.

Rule of Thumb. For the same number of divisions, the magnitude of the midpoint erorr is, roughly, half the magnitude of the trapezoid error. If the concavity is constant, the signs of the midpoint and trapezoid errors are opposite.

(See the discussion on page 352 of the book. There is a good theoretical reason why this "rule of thumb" holds.)

Problem 1. Define the function . Explore the validity of the above rule of thumb by using Maple to calculate TRAP(g,a,b,n), MID(g,a,b,n) for several values of a,b,n. Calculate the corresponding "exact" values, the corresponding errors and quotients. Write a sentence or two describing your results.

As you know, the Rule of Thumb above leads to yet another numerical scheme called the Simpson's rule.

>

Speed of Convergence. The Effects of Doubling n on Errors

The main question when studying a numerical method is how fast the approximation provided by the method converges to the exact value. One way to measure the speed of convergence for our numerical integration schemes is to see how the error decreases when we double the number of divisions n. Hence, for each of our rules of numerical integration we would like to answer the question: by what ratio does the error decrease when n is doubled? Below we investigate this question for the trapezoid rule.

As our test example we shall use the same function and the integral as in the previous section. We shall calculate trapezoid approximations, corresponding errors and ratios for many values of n. It will be convenient to use so called lists . A list in Maple is a sequence of things separated by commas and enclosed by square brackets. The order of elements in a list matters. The latter distinguishes a list of elements from a set of elements. Let's define our list of n values and label it L. Observe that the next entry is obtained by doubling the previous one.

> L:=[5,10,20,40,80,160];

For a given list, L we use the syntax L[i] to denote the i-th element of the list. For example

> L[3]; L[4];

>

Now we want to calculate the trapezoid approximations for our function f, a=0, b=1, for all six values n from list L. Instead of typing six separate commands and labeling the results T1, T2, etc, we can calculate all six approximations at once and put them in a list. We shall label the list Ts. We obtain the list of approximations using the following syntax.

> Ts:=[seq(TRAP(f,0,1,L[i]),i=1..6)];

>

Observe how we use the command " seq ", which stands for "sequence", to evaluate TRAP(f,0,1,L[i]) for all elements in the list L from i=1 to i=6. We surround the " seq " command with square brackets so that the result is another list.

We have obtained a list Ts of values of trapezoid approximations corresponding to the six values of n from list L . Now we shall make a list of errors corresponding to values from the list Ts . The list of errors is labeled " ErrorTs ". Recall that our reference exact value of the integral A is exactA .

> ErrorTs:=[seq(exactA-Ts[i],i=1..6)];

>

Finally, we want to set up the corresponding list of ratios of magnitudes of consecutive errors. Recall that these are the errors corresponding to values of n that are doubled at each step.

> RatioErrorTs:=[seq(abs(ErrorTs[i])/abs(ErrorTs[i+1]),i=1..5)];

>

What do these numbers tell us about the effects of doubling n on the error for the trapezoid rule? It seems that doubling n causes the magnitude of the error to decrease by the factor of 4. It is so in general. More numerical experimentation, with different functions on different intervals leads to the following rule of thumb.

>

Rule of Thumb. For the trapezoid rule, doubling the number of divisions n causes the magnitute of the error to decrease by the factor of four.

Problem 2. In this problem, use the same test function and integral as above.

(a) Investigate the effects of doubling n on the error for the right rule RIGHT(f,a,b,n). Calculate the corresponding errors and ratios for a list of doubling values of n. Use lists in your work with Maple as in the above example. State your answer as a "Rule of Thumb".

(b) Find the errors for the Simpson's rule SIMP(f,a,b,n) for a list of doubling values of n. Note how small those values are. Do not attempt to calculate ratios of errors. Errors are so small that Maple's default setting of acuracy to ten decimal digits does not allow for a good comparison. There is a way to change the default setting, but we are not going to do it at this point.

Problem 3. For the function and integral studied above, what do you think TRAP(f,0,1,320) is? Give your answer based only on the values we computed above for TRAP(40), TRAP(80), and TRAP(160). Then use Maple to verify your answer. How close were you?

Problem 4. For the midpoint rule, what is the effect on the error of tripling n? Experiment with Maple to justify your answer.

Note: If you choose for your experiments a function whose domain is not the whole real line, make sure that the interval of integration is contained in the domain, as it should be. Otherwise, you may cause the program to freeze solid!

Problem 5 . Suppose you have found TRAP(40) and TRAP(80) for a certain integral. Think of a way of combining these to get an estimate of the integral which is more accurate than either one separately. Note that just averaging them might be an improvement, but maybe not. After all, TRAP(80) is likely to be more accurate, so why spoil it by averaging with TRAP(40). But consider that the error, exact - TRAP(40), is about 4 times as large as the error exact-TRAP(80). But we suppose you don't know the value of exact! You might use the preceding sentence as an equation to be solved for "exact." Try and see. Compare your improved estimate with TRAP(40), TRAP(80) and Maple's exact value for the integral A = calculated above.

>

Error Estimates in Terms of n

It turns out that with more careful analysis it can be shown that errors in the trapezoid, midpoint and Simpsons rules satisfy

(1)

(2)

(3) ,

where are constants which depend on the particular function g and the interval [a,b] but do not depend on n . Unfortunately, the constants depend on features of higher derivatives and are usually not available to us. Observe, however, that whatever those constant are, formulas (1)-(3) above mean, roughly speaking, that, as n gets larger and larger, midpoint and trapezoid approximations converge to the actual value of the integral "as fast as" converges to 0, while the Simpson's approximations converge to the actual value as fast as converges to 0, that is, much faster. Do you see how estimates (1)-(3) explain the results obtained above about the effects of doubling n on the error?

>

Graphical Comparison of the Errors for Different Rules

It is always helpful to represent your data graphically. In this section we shall graph errors corresponding to different numerical integration rules. To change pace and to reinforce our skills of working with lists, let's choose a different function to work with and a different list of values for n. Consider

> g:=x->ln(1+x^2);

> B:=Int(g(x),x=0..2); value(B);

>

Again, the actual exact value of the integral B is not very convenient for numerical work. Hence, we shall take as the "exact" value its ten digit approximation.

> exactB:=evalf(value(B));

Define a new list of values for n and label the new list " NL ".

> NL:=[10,20,30,40,50];

Now let's define the corresponding lists of errors for left, right, midpoint, trapezoid and Simpson's rules.

> ErrorLeft:=[seq((exactB-LEFT(g,0,2,NL[i])),i=1..5)];

> ErrorRight:=[seq((exactB-RIGHT(g,0,2,NL[i])),i=1..5)];

> ErrorTrap:=[seq((exactB-TRAP(g,0,2,NL[i])),i=1..5)];

> ErrorMid:=[seq((exactB-MID(g,0,2,NL[i])),i=1..5)];

> ErrorSimp:=[seq((exactB-SIMP(g,0,2,NL[i])),i=1..5)];

>

Observe that each of the above is a list of errors. To plot a list we need the " listplot " command that is contained in the " plots " package. We load the package.

> with(plots):

Since we have worked with the package before, we end our command with a colon so Maple doesn't display the content of the package. We want to plot some of the above lists of errors in one coordinate system, hence we shall label them lp1,lp2 etc and then use the command "display" to plot groups of them together. Each individual command we end with a colon as we don't want individual plots shown.

> lp1:=listplot(ErrorLeft, color=red):

> lp2:=listplot(ErrorRight, color=blue):

> lp3:=listplot(ErrorTrap, color=black):

> lp4:=listplot(ErrorMid, color=green):

> lp5:=listplot(ErrorSimp, color=magenta):

If we try to display all of the above plots in one coordinate system, some of them would be barely visible or not visible at all. That is because errors corresponding to left and right rules are very large in comparison with errors corresponding to other rules, so the graphs of the other rules would merge with the horizontal axis. Let's display them a few at a time.

> display([lp1,lp2]);

>

The above graph gives errors of the left and the right rules. As you see, the " listplot " command plots consecutive points from the list and joins them by segments. Below we display the errors of the midpoint and the trapezoid rules. Notice how much smaller the errors are. Observe also our rule of thumb at work: midpoint errors are of the opposite sign that trapezoid errors and about half of their magnitude.

> display([lp3,lp4]);

>

Finally, we display the error for Simpson's rule. The values are so small that they wouldn't show up on any of the above graphs.

> lp5;

>

Actually, the errors of Simpson's rule approximation are so small and they approach zero so quickly that Maple switched to scientific notation for values on the vertical axis to display the graph. We shall not go into its meaning at this point. Let's just say that errors for Simpson's rule approach zero too quickly to visualize them!

Problem 6. For the function and the integral , list and plot the errors for the left and right rules, and then for the midpoint and trapezoid rules for n=5, 10, 20, 40, 80.

Homework Problems

Your homework for this worksheet consists of Problems 1-6 above.

Note. For this worksheet you will be better off not using the restart command in your homework worksheet. If you do use the command, you have to copy and re-execute all the procedures defining LEFT, RIGHT, and other rules, as well as redefine the function f(x), reload the package with(plots), etc.

MTH 142 Maple Worksheets written by B. Kaskosz and L. Pakula,Copyright 1999.

Last modified August 1999.

>