Integrals 215 problem 04.mws

1. The Trapezoidal Rule.

A Riemann sum to approximate the integral Int(f(x),x = a .. b) is given by Sum(f(c[i])*Delta*x[i],i = 1 .. n) , where Delta*x[i] represents the length of the interval [ x[i-1], x[i] ], c[i] is an arbitrary point in this interval, and the interval [ a, b ] is divided up into n subintervals by the points a = x[0] < x[1] < ...< x[n] = b . If the lengths Delta*x[i] are all equal, so that Delta*x[i] = (b-a)/n , then we have a formula for the points: x[i] = a+i*(b-a)/n . If we further choose the points c[i] according to some reasonable scheme, then we obtain a formula which can be executed by the computer. For example, if we use left endpoints, then c[i] = x[i-1] , or if we use right endpoints, then c[i] = x[i] . The Trapezoid Rule is an average of these two methods, and usually gives a better approximation to the area than either one of these methods. Let us define both the left sums and the right sums in Maple.

> ls:=n->(b-a)/n*sum(f(a+(i-1)*(b-a)/n),i=1..n);

ls := proc (n) options operator, arrow; (b-a)/n*sum...

> rs:=n->(b-a)/n*sum(f(a+i*(b-a)/n),i=1..n);

rs := proc (n) options operator, arrow; (b-a)/n*sum...

We could define the Trapezoid Rule as just being the average of the left and right sum calculations, but we can also define it directly, according to the formula

> ts:=n->(b-a)/(2*n)*(f(a)+2*sum(f(a+i*(b-a)/n),i=1..n-1)+f(b));

ts := proc (n) options operator, arrow; 1/2*(b-a)/n...

It is not difficult to understand that this formula is based on the idea that if you add the left and right hand sums, each term in the sum is counted twice, except for the terms corresponding to x[0] = a and x[n] = b , which only occur once in the sum. Let us use this rule to compute the estimate for the area Int(x^2,x = 1 .. 2)using n = 4 . First we define the function and the endpoints, and then we can use the functions we have defined in Maple to estimate the area.

> f:=x->x^2;a:=1;b:=2;

f := proc (x) options operator, arrow; x^2 end proc...

a := 1

b := 2

Let us compute the left sum, the right sum, their average, and the trapezoid approximation.

> [ls(4),rs(4),(ls(4)+rs(4))/2,ts(4)];

[63/32, 87/32, 75/32, 75/32]

Notice that the trapezoid rule gives exactly the same answer as the average of the left hand and right hand sums.

Submission:

For the function f(x) = sqrt(1+x^4) do the following.

Evaluate the integral Int(f(x),x = 0 .. 1) in Maple, and get a floating point estimate for the answer.

Use the Trapezoid Rule with n = 1000 to estimate the integral.

Finally, find the error in the trapezoidal approximation.

Submission worksheet:

>

2. The error in the Trapezoid Rule.

When `@@`(D,2)(f) is continuous on [ a, b ], the error E[T] in the approximation of the integral using the Trapezoid Rule can be shown to satisfy the estimate

abs(E[T]) <= (b-a)/12*h^2*M ,

where h = (b-a)/n is the length of a subinterval that the interval [ a, b ] is partitioned into by subdividing into n parts, and M is any upper bound for the value of abs(`@@`(D,2)(f)) on this interval. The main difficulty in using this estimate of the error is finding a value of M . If we use the ideas from calculus, we could find the maximum value of abs(`@@`(D,2)(f)) on this interval, and use this exact value of M to give an estimate of the maximum possible error. However, this depends on knowing the function exactly, and in applications that arise in practice, we usually only know the values of the function at certain points by experiment, so we use a more crude estimate for M , taken deliberately to be larger than necessary in order to account for possible experimental error. Also, even if we know the function exactly, finding the upper bound exactly might be very difficult or impossible, so we rely on a choice of M which is easy to arrive at, but is probably larger than necessary.

For example, in Example 3 on page 398, the integral Int(x*sin(x),x = 0 .. Pi) (which Maple can evaluate exactly) is studied. If we define

> f:=x->x*sin(x);

f := proc (x) options operator, arrow; x*sin(x) end...

then we compute the second derivative

> (D@@2)(f)(x);

2*cos(x)-x*sin(x)

Now we could determine exactly where this function has a maximum or minimum on the interval [ 0, Pi ], but it is easy to see that we can choose M = 2+Pi as an upper bound. Let us define a function that will compute the error for the trapezoid rule, given an interval [ a, b ] and an upper bound M , as a function of the number of steps n .

> err:=n->(b-a)/12*((b-a)/n)^2*M;

err := proc (n) options operator, arrow; 1/12*(b-a)...

So for our function, we have

> M:=2+Pi;a:=0;b:=Pi;

M := 2+Pi

a := 0

b := Pi

Suppose that we choose n = 10 steps, then the estimate of the error is

> evalf(err(10));

.1328513704

If we choose n = 100 steps, the error is less than

> evalf(err(100));

.1328513704e-2

Finally, let us define the function, and use the Trapezoid Rule to estimate the integral.

> f:=x->x*sin(x);evalf(ts(100));

f := proc (x) options operator, arrow; x*sin(x) end...

3.141334264

Let us compare this with the value of the integral which Maple (and later in this course, you) can compute exactly.

> evalf(int(x*sin(x),x=0..Pi));

3.141592654

The difference between the Trapezoid Rule and the calculated value of the integral is

> 3.141592654-3.141334264;

.258390e-3

Notice that this difference is smaller than our estimate of the error above.

Submission:

Consider the function f(x) = sqrt(1+x^4) .

Use Maple to find a formula for the second derivative and simplify this formula.

On the interval [0,1], determine a value of M which is an upper bound for the absolute value of the second derivative.

Use this value of M to determine an error estimate for the Trapezoid Rule for n = 1000 .

Compare this value of the error to the actual error that you calculated in problem 1.

Submission worksheet:

>

3. Simpson's Rule.

If you think about the idea behind the Trapezoid Rule, it consists of approximating the function by line segments going through the points [ x[i-1], f(x[i-1]) ], and [ x[i], f(x[i]) ], and computing the area under this approximating curve. The idea of Simpson's Rule is to use parabolas to approximate the function locally. For any three consecutive points on the curve, we can find a quadratic function that passes through these points, and find the integral of this quadratic approximation exactly. Each triple of points corresponds to two intervals if we chop the interval [ a, b ] into smaller subintervals. Thus to use Simpson's rule, we need to divide the interval into an even number of smaller pieces, so we choose the number n of steps to be an even number. The formula which is used to compute the estimate in Simpson's rule is given on page 400 of the text, and is a bit complicated to derive. It can be shown with some work that the formula below, which applies when n is even, is the result, called Simpson's Rule:

> ss:=n->1/3*(b-a)*(f(a)+f(b)+4*Sum(f(a+(2*i-1)*(b-a)/n),i = 1 .. 1/2*n)+2*Sum(f(a+2*(b-a)*i/n),i = 1 .. 1/2*n-1))/n;

ss := proc (n) options operator, arrow; 1/3*(b-a)*(...

Let us apply this rule to estimate the value of Int(5*x^4,x = 0 .. 2) from example 4 on page 400 in the text. We will use n = 4 .

> f:=x->5*x^4;a:=0;b:=2;

f := proc (x) options operator, arrow; 5*x^4 end pr...

a := 0

b := 2

> ss(4);

40/3+2/3*Sum(5*(i-1/2)^4,i = 1 .. 2)+1/3*Sum(5*i^4,...

> evalf(%);

32.08333333

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

32

> 32.08333333-32;

.8333333e-1

Notice how close this estimate is to the real value, even though we used only a very small value of n .

Submission:

For the function f(x) = sqrt(1+x^4) do the following.

First, evaluate the integral Int(f(x),x = 0 .. 1) in Maple, and get a floating point estimate for the answer.

Next, use Simpson's Rule with n = 10 to estimate the integral.

Find the error in the Simpson's Rule approximation.

Compare the error in Simpson's Rule with

n = 10 to the error in the trapezoid rule with n = 1000 .

Submission worksheet:

>

4. Estimating the error in Simpson's Rule.

When `@@`(D,4)(f) is continuous on [ a, b ], the error E[S] in the approximation of the integral using Simpson's Rule can be shown to satisfy the estimate

abs(E[S]) <= (b-a)/180*h^4*M ,

where h = (b-a)/n is the length of a subinterval that the interval [ a, b ] is partitioned into by subdividing into n parts, and M is any upper bound for the value of abs(`@@`(D,4)(f)) on this interval. Just as in the case of the Trapezoid Rule, we use an upper bound M which is not the smallest value possible, but is one which can be found easily.

For example, suppose that we are considering the example from activity 3, the integral Int(5*x^4,x = 0 .. 2) , which we estimate using Simpson's Rule. Let us first define the error estimate, like we did for the Trapezoid rule case.

> err:=n->(b-a)/180*((b-a)/n)^4*M;

err := proc (n) options operator, arrow; 1/180*(b-a...

Now, let us find an estimate for M .

> f:=x->5*x^4;

f := proc (x) options operator, arrow; 5*x^4 end pr...

> (D@@4)(f)(x);

120

Here, we can use the exact value of M =120 in the error estimate.

> M:=120;a:=0;b:=2;

M := 120

a := 0

b := 2

> evalf(err(10));

.2133333333e-2

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

.2133328e-2

Notice that the error estimate is slightly larger than the actual error.

Submission:

Consider the function f(x) = sqrt(1+x^4) .

Use Maple to find a formula for the fourth derivative and simplify this formula.

On the interval [0,1], determine a value of M which is an upper bound for the absolute value of the second derivative.

Use this value of M to determine an error estimate for Simpson's Rule for n = 10 .

Compare this value of the error to the actual error that you calculated in activity 3.

Submission worksheet:

>

>

>