Auxiliary Procedures (execute, but do not change)

>

ApproxIntTable := proc( f, dom, methods, partitions )   local exact, _exact, K, M, hdr, new_row, rows, _hdr, _return, _type;   _hdr := table([right=R,left=L,midpoint=M,trapezoid=T,simpson=S]);   if nargs>4 and convert(args[5],set) subset {`value`,`error`} then     _type := convert(args[5],set)   else     _type := {`value`}   end if;   _return := NULL;   _exact := Int( f, dom );   hdr := [ n, seq( _hdr[m][n], m=methods ) ]:   if member(`value`,_type) then    rows := hdr;    for K in partitions do       new_row := seq( evalf(ApproximateInt( f, dom, method=M, partition=K, output=value )),                       M=methods );       rows := rows, [ K, new_row ];     end do:     _return := _return, Matrix( [rows] )   end if;   if member(`error`,_type) then     hdr := [ n, seq(E[s], s=hdr[2..-1]) ];     exact:=evalf(_exact);     rows := hdr;     for K in partitions do       new_row := seq( evalf(ApproximateInt( f, dom, method=M, partition=K, output=value )                             -exact),                       M=methods );       rows := rows, [ K, new_row ];     end do:     _return := _return, Matrix( [rows] )   end if;   return _return end proc:

>

Lab Overview

Earlier in the semester you spent some time learning a variety of integration techniques - integration by parts, partial fractions, trigonometric substitutions, .... While there are many more techniques that you could learn, the fact is that there are many integrals that cannot be evalauted in terms of our standard functions. Recall that an integral as simple as

could not be evaluated in terms of polynomials, rational functions, or trigonometric functions. In this case, in order to have a name for this antiderivative we created the natural logarithm function:

  .

Another example was seen in Lab 8. When you worked with the normal distribution, you needed to evaluate integrals such as

  .

A new function can be defined to give this a name, but it is not possible to give every integral a name.

The only truly general way to evaluate all definite integrals is numerically. The simplest numerical integration methods involve adding areas of rectangles. In this lab you will see how areas of rectangles can be used to approximate a definite integral. An examination of the errors for a known integral suggest two improved numerical methods, the Trapezoid Rule and Simpson's (or Parabola) Rule.

Example 1: A Graphical Look at

To begin, complete these steps:

Example 2 -A Tabular Look at  

Once again, consider

>	f := 1/x: a := 1: b := 3: II := Int( f, x=a..b ): II = value(II);

The ApproxIntTable  command (defined at the top of this worksheet) can be used to prepare a side-by-side comparison of the three Reimann sum, or Rectangle, Rules for several different partitions. The three Rectangle Rules are

>	Mlist := [ left, right, midpoint ];

and consider the partitions

>	Plist := [ 2, 4, 8, 16, 32, 64, 128, 256, 512 ]:

The corresponding table of approximations to the integral is

>	ApproxIntTable( f, x=a..b, Mlist, Plist, `value` );

Check that these results agree with the ones found in Example 1.

>

The table of errors - the difference between the approximation and the exact value of the integral (ln(3)) - is obtained simply by changing the final argument in the ApproxIntTable  command to `error` :

>	ApproxIntTable( f, x=a..b, Mlist, Plist, `error` );

The important observations in this table are:

•  the left endpoint approximations all overestimate the exact value,
•  the right endpoint approximations all underestimate the exact value,
•  the midpoint approximations all underestimate the exact value, and
•  the midpoint approximations are much more accurate than either the left or right endpoint approximations.

This final observation is most important. To be more explicit, compare the number of subpartitions required to estimate the integral to two decimal places. For the left or right endpoint approximations,  is required; but, for the midpoint approximation,  suffices.

>

>	ApproxIntTable( f, x=a..b, Mlist, Plist, {`value`,`error`} );

>

Example 3 - Trapezoidal Rule

Another observation that can be made from the table of errors in Example 2 is that the errors from the left and right endpoint approximations are very close to being negatives of one another. That is,

  .

A simple average of these two approximations should produce an approximaiton with a much smaller error. The approximation obtained in this manner is called the Trapezoid Rule.

The reason for this name is that the average of the areas of the left and right endpoint rectangles is the area of the trapezoid formed by connecting the left and right endpoints with a straight line. This fact is illustrated in the following figure for

>

N := 2:

subpartitions. (You can change this number, but if it becomes too large you will not be able to distinguish the rectangles and trapezoids.)

>

Pleft := ApproximateInt( f, x=a..b, method=left, partition=N, output=plot,                          title="", boxoptions=[color=blue], showarea=false,                          functionoptions=[color=cyan,style=point] ): Pright := ApproximateInt( f, x=a..b, method=right, partition=N, output=plot,                          title="", boxoptions=[color=green], showfunction=false,                          showarea=false ): Ptrap := ApproximateInt( f, x=a..b, method=trapezoid, partition=N, output=plot,                          title="", boxoptions=[color=red], showfunction=false,                          showarea=false ): display([Pleft,Pright,Ptrap]);

>

To further analyze the Trapezoid Rule, prepare a new table of errors that inludes the errors for the Trapezoid Rule.

>	Mlist2 := [op(Mlist),trapezoid];

>	ApproxIntTable( f, x=a..b, Mlist2, Plist, `error` );

This shows the increased accuracy of the Trapezoid Rule over either the Left or Right Endpoint Rule. But, the Midpoint Rule is always a better approximation than the Trapezoid Rule.

>

Example 4 - Simpson's Rule

Closer inspection of the table of errors in Example 3 reveals the errors for the Midpoint and Trapezoid Rules appears to nearly obey a simple relationship:

  .

If this is true, then a weighted average of  and  should give an improved approximation. The appropriate weighted average is

   =    +   .

The approximation obtained from this formula is called Simpson's Rule. Because this particular weighting amounts to finding the area under a parabola, another name for this method is the Parabola Rule. The following figure adds this parabola to the rectangles and trapezoid. The improved approximation is quite amazing.

>	Psimp := ApproximateInt( f, x=a..b, method=simpson, partition=N, output=plot,                          title="", boxoptions=[color=brown], showfunction=false,                          showarea=false ): display([Pleft,Pright,Ptrap,Psimp]);

Wow! The parabola is almost indistinguishable from the graph of the integrand!! Let's see the errors for this method

>	Mlist3 := [op(Mlist2),simpson];

>	ApproxIntTable( f, x=a..b, Mlist3, Plist, `error` );

Notice that Simpson's Rule with  provides an approximation that is accurate to two decimal places. (Compared with  for the Midpoint and Trapezoid Rules and  for the Left and Right Endpoint Rules.)

Looking towards the bottom of the column of errors for Simpson's Rule, these errors are not really zero.  Maple reports 0.  because the actual error is smaller than the smallest floating-point number with 10 decimal digits, i.e., . To force Maple to work with smaller floating point numbers, change the Digits  parameter to the number of digits Maple should keep, say 14,, then re-create the table of errors.

>	Digits := 14: ApproxIntTable( f, x=a..b, Mlist3, Plist, `error` ); Digits := 10:

To conclude, note that each time the number of subintervals is doubled, Simpson's Rule provides one additional decimal digit in the approximation. (See Lab Question 2.)

>

Lab Questions