Math 166
 Computer Lab # 4
 Numerical Integration

Goal:  To approximate definite integrals using left and right   Riemann sums, the midpoint rule, and the trapezoidal rule; to approximate integrals when the integrand is given discretely.

New Mathematica  (built-in) commands introduced:  Sum 

New commands we defined:  leftRiemannsum, rightRiemannsum, midpointsum, trapezoidalsum, plottrapezoid

Introduction

In your study of calculus you have learned to evaluate certain  integrals.  However, for many integrals it is not easy, or just not possible, to find an elementary antiderivative of the integrand so that we may evaluate the integral. In other instances the function representing the integrand is not given explicitly but rather by a set of discrete data.  The purpose of this lab   is to look at some techniques for approximating the value of a definite integral when the integrand is given either explicitly or as a set of discrete data.  Note that Mathematica  has a built-in command 
             NIntegrate[f,{x,xmin,xmax}] 
for approximating the value of an integral but it does not allow us to choose the particular method it uses.  We want to approximate the values of definite integrals by methods of our choice.

Right and Left Riemann Sums

When first studying definite integrals the concept of Riemann sums is used to approximate definite integrals.  This was a somewhat tedious task which involved repeated sums of areas of rectangles.  This process can be easily automated with Mathematica  using the Sum command which was introduced in class.  

If we have an explicit function given for the integrand and we are using equally spaced intervals,  then we can define a function in Mathematica  to automatically calculate the left or right Riemann sum for a given function, specified endpoints, and a given number of intervals.  We have the following definition for a left Riemann sum which outputs the answer in the variable approx.

Clear[leftRiemannsum]
leftRiemannsum[f_,{x_,a_,b_,n_}]:=
 Block[{delx,i},
   delx=(b-a)/n;
   approx=N[ Sum[(f /.x->a+delx*i)*delx,{i,0,n-1}] ] ]

Note that the arguments to this function are similar to the command Integrate except we must now include the number of equal subintervals that we are using.  Specifically, the arguments are (i) f , the integrand, (ii) x, the independent variable, (iii) a, the lower endpoint, (iv) b, the upper endpoint, and (v) n, the number of equal subintervals.  The  effect of the Block command  is to make the variables in the list local variables.  In our case it is the variables delx and i. The length, delx,  of each of the equal subintervals is calculated in the usual way of dividing the total length, b-a, by the number of subintervals.  The notation /. in the definition of approx is a replacement operator so that the notation  expr /. lhs -> rhs  replaces the term lhs with rhs in expr; that is, f is evaluated at x=a+delx*i.  We can define the function for the right Riemann sum analogously.

Clear[rightRiemannsum]
rightRiemannsum[f_,{x_,a_,b_,n_}]:=
 Block[{delx,i},
   delx=(b-a)/n;
   approx=N[ Sum[(f /.x->a+delx*i)*delx,{i,1,n}] ] ]

Example 1.  Approximate the integral  Ü13 e-x^2 dx using a left and right Riemann sum with n=20,40,60,80 equal subintervals.  Output n and  your approximations into table format with headings.  Explain the results.

To simplify the calculations one can use the Table command with its first argument a list including the number of intervals, the value of the left Riemann sum and the value of the right Riemann sum for each number of intervals.

Clear[f,x]
f[x_]:= E^(-x^2)
Table[ {n,leftRiemannsum[f[x],{x,1,3,n}],
          rightRiemannsum[f[x],{x,1,3 ,n}]}, {n,20,80,20}];
TableForm[%,TableHeadings->{None,{"n","lhs","rhs"}}]

As the number of subintervals is increased, the number of significant digits in the approximation is increasing.  For example, between n=20 and n=40 the approximation is changing in the second digit; between n=40 and n=60 it is changing in the third digit, etc.  Also note that since the function e(-x^2) is decreasing on the interval [1,3] we expect the right Riemann sum to be an underestimate for the area and the left Riemann sum an overestimate as indicated by our data.

Note that in the above example this is an integrand for which an antiderivative can not be easily found.  However, we can use the NIntegrate command to see how our approximation compares with Mathematica's.  In this case Mathematica is using a more accurate method so we expect a better answer.

NIntegrate[E^(-x^2),{x,1,3}]

Exercise 1.  The midpoint rule uses rectangles whose height is determined by the value of the function at the midpoint of the subinterval.  Modify the function for the right or left Riemann sum to define a function which computes the midpoint rule for n equal subintervals.  Use the midpoint rule to approximate the integral in Example 1 using n=20,40,60,80 equal subintervals.

Note that we have simply evaluated the function at the midpoint of the interval, i.e., at xi+delx/2

Clear[midpointsum]
midpointsum[f_,{x_,a_,b_,n_}]:=
 Block[{delx,i},
   delx=(b-a)/n;
   approx=N[ Sum[(f /.x->a+delx*i+delx/2)*
            delx,{i,0,n-1}] ] ]

Clear[f,x]
f[x_]:= E^(-x^2)
Table[ {n,midpointsum[f[x],{x,1,3,n}] }, {n,20,80,20}];
TableForm[%,TableHeadings->{None,{"n","midpt" }}]

Trapezoidal Rule

We have seen in class that we can use a trapezoid instead of a rectangle to approximate the area under a curve.   This is shown graphically below.

Clear[plotf]
  plotf=Plot[x^4+1,{x,0,2.5},PlotStyle->
      RGBColor[1.000, 0.144, 0.077],
      DisplayFunction->Identity,
      Ticks->None,PlotRange->All,
      PlotLabel->"Trapezoidal Rule" ]
  trap= Line[{{1,0},{1,2},{2,17},{2,0}}];
  Show[plotf,Graphics[trap],Graphics[Text["a",{1,-2}]],
  Graphics[Text["b",{2,-2}]],
  Graphics[Text["f(x)",{2.3,22}]],
     DisplayFunction->$DisplayFunction ]

Recall that the trapezoidal rule can be written in either of the following ways.

Trapezoidal Rule:
   Üab f(x) dx â ( [f(x0)+2f(x1)+2f(x2)+...+2f(xn-1)+f(xn)]*êx/2  
where êx=(b-a)/n;
or
   Üab f(x) dx â average of left and right  Riemann sums.

Thus we can easily define a function for the trapezoidal approximation in terms of the right and left Riemann sums.

Clear[trapezoidalsum]
trapezoidalsum[f_,{x_,a_,b_,n_}]:=
  N[ ( leftRiemannsum[f,{x,a,b,n}]+
      rightRiemannsum[f,{x,a,b,n}] ) /2 ]

Exercise 2.  Use the trapezoidal rule to approximate the integral  Ü13 e-x^2 dx using  n=20,40,60,80 equal subintervals.  How does this compare with the values you get using the right and left Riemann sums with the same number of intervals?  How does this compare with the results from the midpoint rule?

Clear[f]
f[x_]:=E^(-x^2)
Table[{n,trapezoidalsum[f[x],{x,1,3,n}] },{n,20,80,20}];
TableForm[%,TableHeadings->{None,{"n","trap" }}]

The error in the trapezoidal rule is clearly much smaller than either the right  or left  Riemann sums and is comparable with  the error obtained using the midpoint rule.

Plotting our Approximations

It's always informative to see our approximations graphically as well as in a table of data.  The following command allows you to plot a function f(x) from x=a to x=b using n trapezoids with base (b-a)/n so we can see graphically the effect of adding more and more trapezoids.    The form of the command is the same as for the approximations:      plottrapezoid[ f[x], {x,a,b,n} ]   .

Clear[plottrapezoid]
plottrapezoid[f_,{x_,a_,b_,n_}]:=
  Block[ {delx,xi,trap,plotf,heights,approx},
  plotf=Plot[f,{x,a,b},PlotStyle->
      RGBColor[1.000, 0.144, 0.077],
      DisplayFunction->Identity];
  delx=N[(b-a)/n];
  xi[i_]:=a+delx*i;
  heights=Table[f/.x->xi[j],{j,0,n}];
  trap=Table[ Line[{ {xi[j],0 },{xi[j],heights[[j+1]]},
   {xi[j+1],heights[[j+2]]}, {xi[j+1],0}}],{j,0,n-1}];
  Show[plotf,Graphics[trap],
     DisplayFunction->$DisplayFunction,
     PlotLabel->   "trapezoidal rule" ]]

Exercise 3.   Consider the following function.

Clear[g]
g[x_]:=400*x^5-900*x^4+675*x^3-200*x^2+25*x+0.2

We want to approximate the integral Ü0.8 g(x) dx using n=2,4,8,16,32 equally spaced trapezoids.  Make a table of n, your approximation and the exact error (that is, use the command Integrate to exactly evaluate the integral  and give the absolute error).

exact=Integrate[g[x],{x,0,.8}] 
Table[{n,trapezoidalsum[g[x],{x,0,.8,n}],
  Abs[exact- trapezoidalsum[g[x],{x,0,.8,n}]] },
  {n,2,32,2}];
TableForm[%,TableHeadings->{None,{"n","trap","error"}}]

Exercise 4.   For the same function as in Exercise 3, use the command plottrapezoid to plot your approximations.
Animate the four plots using the Animate Selected Graphics from the Graph menu.  You may want to reset the number of frames per second. This is done by going to the menu Graph and selecting Animation; a window will pop up and you can reset the frames per second -- I suggest to 1 per second.

plottrapezoid[g[x],{x,0,.8,2}]
plottrapezoid[g[x],{x,0,.8,4}]
plottrapezoid[g[x],{x,0,.8,8}]
plottrapezoid[g[x],{x,0,.8,16}]
plottrapezoid[g[x],{x,0,.8,32}]

Handling Problems with Discrete Data

If the function is not given explicitly, we can not use the Mathematica  functions such as leftRiemannsum or trapezoidalsum that we have defined to approximate the definite integral.  Instead, we must return to using the Sum command or alternately define a new function.  If the integrand is represented as a set of discrete data we handle this in Mathematica  by putting the data in lists.  If the function is given graphically, then we must obtain the data from the graphs and input it into lists.  The following example illustrates a problem where we are given a set of discrete data representing the integrand.

Example 2.  A car accelerates from 0 to 60 mph in 12 seconds.  Its velocity is recorded during this 12 second period and is given in the table below.  Use the trapezoidal rule to estimate the distance in feet the car traveled in the 12 seconds.

t  (in sec)     1     3     4     5      7     9     11     12
v(t) (in mph)   4    14    21    27     39    47     57     60

We want to approximate the definite integral of v(t) from t=0 to t=12.  However, our integrand is given by the above set of discrete data instead of by an explicit formula.  Note that the data is not given at equally spaced time intervals.  We enter our data into 2 lists called time & vel with 8 elements in each list. (Recall that lists in Mathematica are enclosed in 
{ } and to access an element in a list we use [[ ]]; e.g., the fourth element in the list time is obtained by typing time[[4]].) Since the time is not given at equally spaced intervals, we must multiply the height of each rectangle (given by the velocity) by the change in time; e.g., for the first interval the change in time is (3-1)=2 seconds.  Note also that the time is given in seconds and the velocity is given in mph so that we must convert the velocity to feet per second to obtain our answer in feet; i.e., multiply velocity by 5280ft/3600sec.

time={1,3,4,5,7,9,11,12};
vel={4,14,21,27,39,47,57,60};
ftps=5280/3600;
lhs=N[ Sum[ftps*vel[[i]]*(time[[i+1]]-time[[i]]),{i,1,7}] ]
rhs=N[ Sum[ftps*vel[[i+1]]*(time[[i+1]]-time[[i]]),{i,1,7}]  ]
distance=(lhs+rhs)/2

Thus the approximation to the distance traveled using the trapezoidal rule is 548.533 feet.

Problems to be handed in

1. A function  called the integral logarithm  is defined by
      Ü2 x g(t) dt    for xÑ 2   where g(t) = 1/ln t.
This function occurs in analytic number theory where it is proved that it is a good approximation to the number of primes less than or equal to x.  However, we do not know an antiderivative of g(t) so we approximate the integral using a Riemann sum.

a.  Use the midpoint rule with n=100, 200, 300, 400 to approximate the number of primes æ 5000.

b.  The built-in command Prime in Mathematica  gives the ith prime number; for example, Prime[2] gives the second prime number listed in ascending order, i.e., 3.  The Table command can be used to print out a list of  prime numbers. 
Use the Prime command to find the exact number of prime numbers æ5000 and compare with your approximation in (a). ( Think!  You don't want to type Prime[5000] since it will list the first 5000 prime numbers.) 

c.  Repeat (a) and (b) to approximate the number of primes
æ10,000.

2.  We want to approximate  Ü13 ( 1 / (1+x2) ) dx for n=20,30,40,...80  equal subintervals using   (i) a right Riemann
sum, (ii) a left Riemann sum, (iii) the midpoint rule, and (iv) the trapezoidal rule.

a.  Make a table of the value of n and each of the approximations (i)--(iv).  Label each column using the TableHeadings option of the TableForm command.  Plot your approximations using the trapezoidal rule.

b.  Find the exact value of the integral using the command Integrate. Make a table of the value of n and the error in each of the approximations (i)--(iv).  Label each column.  Discuss your results.  In particular, do you see any pattern in the error as n is doubled, i.e., the length of the subinterval is halved?

c.  Plot the function 1/(1+x2)  from x=1 to x=3.  State whether the function is increasing or decreasing and discuss its concavity.  Use this to determine if each approximation in (a) is an underestimate or an overestimate.

3. Suppose that we are given the air temperature every three hours for a particular day in September for Ames and want to use this data to estimate the average temperature over the 24 hour span.  (Recall from page 284 in your text that the average value of a function f defined on an interval [a,b] is given by
      (1/(b-a) Üab f(x) dx . )
The following lists give the hours of the day and the corresponding temperature in Fahrenheit degrees; in addition, we have defined another list which gives a list of coordinate pairs of the form (hour,temperature) for plotting purposes.  Use this data to obtain an estimate for the average temperature using the trapezoidal rule.  Plot the data.

hour={0,3,6,9,12,15,18,21,24};
temp={60,56,53,65,74,76,71,67,60};
data=Table[ {hour[[i]],temp[[i]]},{i,1,9}];

4. A sketch of a fairway on a golf course is given in the figure below; all measurements are given in feet.  Use the trapezoidal rule with n=10 equal subintervals to estimate the amount of fertilizer needed to fertilize this fairway if one pound of fertilizer covers 500 square feet.

h=100;
data={ {0,h},{0,h+5},{100,h+40},{200,h+30},{300,h+25},
{400,h+15} ,{500,h+12},{600,h+16},{700,h+34},
{800,h+40},{900,h+50},{1000,h+10},{1000,h},
{900,30},{800,35},{700,34},{600,16},{500,7},
{400,5},{300,30},{200,45},{100,80},{0,h} };
ListPlot[data,PlotJoined->True,GridLines->
{{100,200,300,400,500,600,700,800,900,1000},Automatic},
AxesLabel->{" feet","feet"}]