The Trapezoidal Rule

Copyright 2000

Department of Mathematics

University of Georgia

Athens, Georgia

Carol W. Penney

Robert Rumely

John Gosselin

The definite integral int(f(x),x = a .. b)  is defined to be the limit of Riemann sums as n approaches infinity .  Because of the difficulty of finding the limit of the Riemann sum we usually do not attempt to evaluate an integral directly from this definition.  We generally prefer to use the Fundamental Theorem of Calculus, which enables us to find the exact value of     int(f(x),x = a .. b)   by calculating F(b)-F(a), where F is an antiderivative of  f  on [a,b].  Occasions exist, however, in which we prefer to approximate the value of an integral numerically rather than using the Fundamental Theorem to evaluate it exactly.  Two such occasions are:

 

It is amazing that there are methods enabling us to accurately approximate integrals with a relatively small amount of computation.   In this project we study one of these methods.  Our focus will be on understanding why they work and how we can implement them.

Consider the example

   Int(x^2*exp(-x),x = 0 .. 5)

which is an integral that can actually be computed exactly.  We will approximate it by both Riemann sums and a new type of sum that involves trapezoids - hence the name "Trapezoidal Rule".

 

To approximate an integral    Int(f(x),x = a .. b) ,   these methods work by dividing the interval  [a,b]  into  n  equal pieces  and adding up an appropriate combination of values of   f(x)  at the partition points.   Naturally the larger  n  is, the more accurate the result should be (subject to round-off error, of course).  Offsetting this desire for accuracy is the fact that larger values of  n  require more calculations.  We start off with the case  n = 8. We first consider Riemann sums. We will only considerright Riemann sums.

The Riemann Sum Rule

The " Riemann Sum Rule " for approximating   Int(f(x),x = a .. b)    is a procedure you have already encountered:  subdivide  [a,b]  into  n  equal parts (the regular partition), and compute the Riemann sum corresponding to a selection, here given by the right-hand endpoints.  Since the  integral is defined to be a limit of Riemann sums, for  n   sufficiently large the value of the Riemann sum should be reasonably close to that of the integral.

 

Actually we do not advocate using this "Rule" in practice:  the other methods we will introduce later are better. You may be wondering, why not use the midpoint Riemann sum, rather than the right endpoint one? The "Midpoint Rule" obtained by selecting the midpoint of each subinterval is better than the "Riemann Sum Rule", but it has other disadvantages which will become clearer later.  For the moment, please bear with us, and view the "Riemann Sum Rule" using the selection of the right-hand endpoints as an example, setting a pattern for our investigation of the other rules.

Graphical Illustration of the Riemann Sum Rule

>    restart:

Define the function, interval, and  n:

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

>    a:=0.0;

>    b:=5.0;

>    n:=8;

The width of each of the 8 subintervals will be called  h.  The value of  h  is as  follows:

>    h:=(b-a)/n;

For i=0, 1, 2, . . . ,n, define x(i):

>    x:=i->a+i*h;

Make the selection of   xstar(i)  as the right endpoint.

>    xstar:=i->x(i);

Define the ith rectangle by using the rectangle command, which defines the rectangle by listing a pair of opposite vertices:

>    riemannrectangle:=i->rectangle([x(i-1),0],[x(i),f(xstar(i))],  color=red);

Let rectangles  be the list of rectangles approximating the area of the region under the graph of y=f(x).

>    rectangles:=[seq(riemannrectangle(i),i=1..8)]:

Load plottools , and plot the graph of  y=f(x)  and the Riemann rectangles.

>    with(plots):with(plottools):

>    display(rectangles,plot(f(x),x=0.0..5.0,y=0.0..0.6,
color=navy,thickness=2,title=`Riemann Sum Approximation`));

>   

The Formula for the Riemann Sum Rule

Now we calculate the value of the Riemann sum:  the i^th  rectang %? le has width  h=(b-a)/n   and height    y[i] = f(x(i)) .  We usually write the right Riemann sum in the form sum(f(x(i))*h,i = 1 .. 8)  .  It is convenient for this project to rewrite the sum in the form h*( y[1]+y[2] + ... + y[8] ). In general if we were using n  subintervals, the right Riemann sum could be written in the form h*( y[1]+y[2] + ... + y[n] ).

 

This is an important formula to remember.  Now we will evaluate this sum (the area of the red rectangles in the picture above).  First define  y(i) to be equal to  f(x(i)):  

>    y:=i->f(x(i));

Now write the sum given above:

>    RRApprox:=h*sum(y(i),i=1..8);

The Error in the Riemann Sum Approximation

Let us compare this approximation to the value obtained with the integrate  command:

>    TrueValue:=evalf(int(f(x),x=a..b));

Error:  The error  in this approximation is defined to be the absolute value of the difference between the value of the integral and the value of the approximation. This error is:

>    RiemannError:=evalf(abs(RRApprox-TrueValue));

To calculate the percentage error, divide the error by the actual value and multiply by 100:

>    100*(RiemannError)/TrueValue,`percent`;

>   

>   

In this example we used only 8 rectangles, yet the approximation was within  3%  of the true value.   Does this surprise you?  

The Trapezoidal Rule

If you look at the picture of the Riemann approximation, you will see that the essence of the Right-Hand Riemann Sum Rule is that we approximate y=f(x), on each subinterval, by a horizontal line--the graph of a constant function.  It would seem reasonable to expect that a straight (not necessarily horizontal) line running between two successive points on the graph, that is from point  (x(i-1),f(x(i-1)))  to point  (x(i),f(x(i))),  would be a better approximation to the function, since it approximates the slope as well as the height of the function. That is, an improvement method of approximating f(x) on each subinterval by a constant function would be to approximate f(x) on each subinterval by a general linear function.  This idea results in trapezoidal areas rather than rectangular areas   approximating the area under the graph, and hence the name Trapezoidal Rule.  

Graphical Illustration of the Trapezoidal Rule

>    restart:

Define the function f(x)  and the original values of  a, b, and n:

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

>    a:=0.0;

>    b:=5.0;

>    n:=8;

The width of each trapezoid is h:

>    h:=(b-a)/n;

The x(i) and y(i):

>    x:=i->a+i*h;

>    y:=i->f(x(i));

Now load the appropriate packages:

>    with(plottools):with(plots):

Define the ith trapezoid as a polygon, by listing its vertices:

>    trapezoid:=i->polygon(
[[x(i-1),0],[x(i-1),y(i-1)],[x(i),y(i)],[x(i),0]],
color=aquamarine):

Constructing the sequence of trapezoids:

>    trapezoids:=[seq(trapezoid(i),i=1..8)]:

>    display(trapezoids,plot(f(x),x=0.0..5.0,y=0.0..0.6,
color=navy,thickness=2, title=`Trapezoidal Approximation`));

>   

The area of the trapezoids appears to give a more accurate approximation to the definite integral than the Riemann sum approximation did.  To determine whether this idea is useful, we need to find out whether calculating the area inside these trapezoids is significantly more difficult than calculating the Riemann sum.

The Formula for the Trapezoidal Rule

Let's evaluate this trapezoidal sum. Consider the following trapezoid, which is the second trapezoid of our approximation.  Let it have width  h,  left height  y(1),  and right height y(2) .  It consists of a rectangle surmounted by a triangle, so its area is the sum of the area of the pink rectangle and the aquamarine triangle displayed below.

[Maple Plot]

The area of the rectangle is   y[1] *h;  the area of the triangle is (1/2)*( y[2]-y[1] )*h.  We first restart to remove the definition of h:

>    restart:

The area of the trapezoid is therefore

>    y1*h+(1/2)*(y2-y1)*h;

Simplify and factor, to obtain the formula for the area of our second trapezoid:

>    factor(simplify(y1*h+(1/2)*(y2-y1)*h));

This calculation would be the same for each trapezoid; adding the areas of our 8 trapezoids gives:

>    factor(sum((1/2)*h*(y[i-1]+y[i]),i=1..8));

which is the formula for the Trapezoidal Rule approximation, for n=8.   Notice that it has only one more term than the Right-Hand Rule sum with n=8, so is not significantly more difficult to compute.   The sum in parentheses is called a weighted sum ;  the heights of the function, for the trapezoidal approximation, are weighted-- the weights c[i]  being the coefficients 1,  2,  2,  2,  2, . . ., 2,  1, where the number   y[i] = f(x(i))   is the height at the   i^th   point.  

 The formula for the Trapezoidal Rule approximation    A[Trap](n)    to    Int(f(x),x = a .. b)    is therefore

A[Trap](n)    =   (1/2)* h*( Sum(c[i]*f(x(i)),i = 0 .. n)  )   =   (h/2)*( y[0]+2*y[1]+2*y[2] + ... + 2*y[n-1]+y[n] )

This is an important formula to remember.  Now we will evaluate this sum for our example (the area of the aquamarine trapezoids in the picture above).  We first remember to define everything again:

>    restart:

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

>    a:=0.0;

>    b:=5.0;

>    n:=8;

>    h:=(b-a)/n;

The x(i) and y(i):

>    x:=i->a+i*h;

>   

>    y:=i->f(x(i));

Notice that in the following command multiply h/2 by this sum:   we first add y(0), then twice the sum of y(1) through y(7), and finally add y(8), to get our weights 1, 2, 2, ,2, 2, 2, 2, 2, 1:

>    TrapApprox:=(h/2)*(y(0)+2*sum(y(i),i=1..7)+y(8));

>   

The Error in the Trapezoidal Approximation

Again let us compare this approximation to the value obtained with the integrate command:

>    TrueValue:=evalf(int(f(x),x=a..b));

  • The error of this approximation is:

>    TrapError:=abs(TrapApprox-TrueValue);

For a percentage error of:

>    100*(abs(TrapApprox-TrueValue))/TrueValue,`percent`;

>   

How does this error compare with the Riemann sum error from the previous section, using the same function and the same number, 8,  of subintervals?

Project

In your project you are to compare the trapezoidal approximations and right Riemann approximations with the aid of a spreadsheet. Let m  denote the number of letters in your first name. Define the function f(x) = 1.0/(1+m*x^3)  on the interval [1,2]. Define the the right Riemann approximations and trapezoidal approximations in terms of n. also define the corresponding errors as functions of n . You should be able to set up these formulas from the work above. Remember that when the interval [a,b] is divided into n subintervals, the Riemann sum involves n  values of  f whereas the trapezoidal rule involves n+1  values of f .

>    restart:

>    f:=x->;

>    a:=;

>    b:=;

>    TrueValue:=;

>    deltax:=n->;

>    xstar:=(i,n)->;

>    RRApprox:=n->;

>    RRError:=n->;

>    TrapApprox:=n->;

>    TrapError:=n->;

Insert a spreadsheet into your report similar to the one below. Check that the number of decimals being displayed has been set to 8 with 10 digits being used for calculations. To check this, right click in one of the blue cells and select properties. Make the necessary changes and then click OK.

n RRApprox TrapApprox RRError TrapError
20



40



60



80



100



NULL



Include the following in your report

The Most Common Maple Commands

Academic Honesty Statement:

Place the following statement (by copying and pasting) at the end of your report and sign it in ink.  Your instructor will not grade your report unless this signed statement appears at the end of your report.

I understand that I may work with others if I give them credit in this statement.  I also understand that I am required to write my report--that to copy all or part of someone else's report or to allow someone else to copy all or part of my report constitutes plagiarism, which is a serious violation of academic honesty.

I worked with (replace this parenthetical remark with first and last names of those with whom you worked)   on this project.  I wrote my own report.  I did not copy any of this report from anyone else and I did not allow anyone else to copy any of this report.

Signed: