135lab2.htm

Math 135 Calculus I - Lab 5 - Numerical Integration Methods I

Math 135 - Calculus I
University of Colorado at Colorado Springs
Lab 5: Numerical Integration Methods I

We have seen two basic approaches for determining the value of a definite integral: the left or right hand Riemann sums. While the underlying ideas for each approach are the same, the obvious distinction between the methods gives us a means by which we can estimate integrals to an arbitrary level of precision - in most cases. Complications may arise when we try to "trap" the exact value between the left and right sums. In this lab you will be asked to investigate some of these complications. Our main tools will be Maple routines that are contained in the "student" package. Since these routines are contained in a package, Maple does not automatically load them when a session is started. You have to instruct Maple to do so:

> with(student);

Suppose we wish to estimate the value of ò ₂⁴x²dx using a left-hand Riemann sum with ten subdivisions. Since is an increasing function on [2,4], we would expect the left-hand sum to give us an underestimate for the true value of the integral. This observation is readily made if we use "leftbox" to generate a picture of the Riemann sum:

> leftbox(x^2,x=2..4,10);

When using "leftbox" command (or any of the other three routines) we first specify the function, "x^2", then the interval, x=2..4, and finally the number of subdivisions. The sum is computed with "leftsum":

> leftsum(x^2,x=2..4,10);

Notice that Maple does not return a numerical value. It gives us a symbolic, or exact, value for the sum. We must tell Maple to give us a floating point (or decimal) answer by using the "evalf" command:

> evalf(%);

The percent sign in the evalf command is Maple shorthand that tells the computer to use the result of the preceding command. Rather than perform two separate commands, we can combine this into one command:

> evalf(leftsum(x^2,x=2..4,10));

Provided that the integrand is monotonic (increasing or decreasing) over the interval in question, we know that the exact value must lie somewhere between the left and right sums. If we increase the number of subdivisions to the point that the difference between these sums (in absolute value) is smaller than the amount of allowable error, we are assured that the error in our approximation is no greater than this difference. But suppose that the integrand is both increasing and decreasing over the interval in question, how do get an accurate estimate in this case?
Suppose we wish to determine the value of with an error of less than .001, i.e., we want to be accurate to at least three decimal places. We could try to compare the left and right sums using, as an initial guess, say, five subdivisions:

> evalf(leftsum(exp(-x^2/2),x=-2..2,5));

> evalf(rightsum(exp(-x^2/2),x=-2..2,5));

Notice that the left and right sums yield identical values. Does this mean that we have found the exact value of the integral using only five subdivisions? Definitely not. This is obvious if we look at the sums with "leftbox" and "rightbox":

> leftbox(exp(-x^2/2),x=-2..2,5);

> rightbox(exp(-x^2/2),x=-2..2,5);

It turns out that for our example, the left and right sums will always yield identical values for number of subdivisions since the function and the interval over which we are integrating are symmetric with respect to the y-axis. However, we have no means of knowing how closely that value is to the 'true' value of the integral. It is not to difficult to work around this problem however. First note that our function is increasing on the interval [-2,0] and decreasing over [0,2]. If we break the integral into two integrals, , we can then trap the value of each integral between left and right sums. Why? By assuring that our estimate for each integral has an error of less than .0001, then the total error cannot be greater than .0002 and our estimate for will be accurate to at least within .001.

Knowing that the integrand is monotonic over the interval of integration means that we can trap the actual value of the integral between the left and right sums. At this point it would be possible to estimate the value to any required accuracy by trial and error. We could simply compare the left and right sums for various numbers of subdivisions until we found an 'n' where the difference between the sums was small enough. However, it might be a more efficient use of our time to find the number of subdivisions necessary before we start asking Maple to evaluate left and right-hand sums.

Recall that the difference between the left and right Riemann sums (in absolute value) is given by the formula

| left sum - right sum | = where.

Substituting for yields the relation | error | . We can control the error in our approximation by making the quantity on the right small enough. Suppose we want to estimate with an error of less than .0001 as in the previous example. In this case, and b = 0. We want to make the quantity since this would insure that the error would also be less than .0001. First we define our function f(x) in Maple using the following command:

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

Now Maple 'knows' of f as a function of x.

> evalf(f(-2));

Next we substitute the quantities into the expression and solve for n. First, we ask Maple to evaluate {|f(-2)-f(0)|(2)}/n, then we ask it to solve to find when this expression is <.0001.

> evalf(abs(f(-2)-f(0))*2/n);

> solve(% <.0001);

Maple gives us two solutions, the numbers less than 0 and greater than 17293.29434. Since we can't have a negative number of subdivisons, n must be greater than or equal to 17294. Notice that since the difference between the left and right sums is less than .0001 for this many subdivisions, we need only calculate one of the sums and we are guaranteed it will approximate the integral within .0001.

> evalf(leftsum(f(x),x=-2..0,17294));

Exercises
1. Compute the following integrals to two decimal places. For each integral, determine the minimum number of subdivisions required to achieve the desired level of accuracy.

a)	c)
b)	d)

2. Using the ideas detailed above, approximate the value of to at least three decimal places. It might be useful to look at a plot of the function (Remember the plot command: "plot(sin(x),x=0..3*Pi/2);").