Riemann, Trapezoids, and Simpson

Off[NIntegrate::"ploss"];Off[NIntegrate::"ncvb"];Off[$MaxExtraPrecision::"meprec"];Clear[riemannsum];
riemannsum[fnc_,tind_,to_,tf_,nint_,rtorlt_,fncicin_:Automatic]:=Block[{},If[fncicin===Automatic,fncic=0,fncic=fncicin];f[t_]=fnc/.tind->t;F[t_]=NDSolve[{y'[t]==f[t],y[to]==fncic},y[t],{t,to-0.1` (tf-to),tf+0.1` (tf-to)}][[1,1,2]];h=((tf-to)/nint);If[rtorlt=="left",discvalues=Table[N[f[t]],{t,to,tf-h,h}],discvalues=Table[N[f[t]],{t,to+h,tf,h}]];Fdiscvalues=Table[{to+k h,fncic+(∑_(i=1)^k)discvalues[[i]] h},{k,1,nint}];Fdiscvalues=Prepend[Fdiscvalues,{to,fncic}];fvalues=Table[f[t],{t,to,tf,((tf-to)/500)}];Fvalues=Table[F[t],{t,to,tf,((tf-to)/500)}];Fmax=Max[Flatten[{Fvalues,Table[Fdiscvalues[[i,2]],{i,1,nint}]}]];Fmin=Min[Flatten[{Fvalues,Table[Fdiscvalues[[i,2]],{i,1,nint}]}]];fmaxt=Max[fvalues];fmint=Min[fvalues];fmin=fmint;fmax=fmaxt;changesign={};Do[If[Sign[discvalues[[k]]]≠Sign[discvalues[[k+1]]],changesign=Append[changesign,to+k h]],{k,1,nint-1}];changesign=Prepend[changesign,to];changesign=Append[changesign,tf];size=7;steprat=0.02`;hrange=If[Fmax≤0,{1.1` Fmin,0},If[Fmin≥0,{0,1.1` Fmax},{1.1` Fmin,1.1` Fmax}]];frange=If[fmax≤0,{1.1` fmin,0},If[fmin≥0,{0,1.1` fmax},{1.1` fmin,1.1` fmax}]];p1=Plot[f[t],{t,to,tf},PlotRange->{{to,tf},frange},AspectRatio->1,PlotStyle->{RGBColor[0,1,0],Thickness[0.008`]},ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity];p2t={};Do[If[discvalues[[k]]>0,fillcolor=RGBColor[1,0,0],fillcolor=RGBColor[0,0,1]];pinter=Plot[discvalues[[k]],{t,to+(k-1) h,to+k h},PlotRange->{{to,tf},frange},PlotStyle->{RGBColor[0,1,0],Thickness[0.008`]},ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity,Filling->Axis,FillingStyle->fillcolor];p2t=Append[p2t,pinter],{k,1,nint}];p2t=Append[p2t,p1];p2=Show[p2t,ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},AxesLabel->{x,fnc},DisplayFunction->Identity];p3=Plot[F[t],{t,to,tf},PlotRange->hrange,ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity];p4=ListPlot[Fdiscvalues,PlotRange->hrange,PlotStyle->{PointSize[0.02`],RGBColor[0.398444`,0.847669`,0.347662`]},ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity];p5=Show[{p4,p3},ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity];p6t={};Do[If[discvalues[[k]]-f[to+(k-0.5`) h]>0,fillcolor=RGBColor[1,0,0],fillcolor=RGBColor[0,0,1]];pinter=Plot[{discvalues[[k]],f[t]},{t,to+(k-1) h,to+k h},PlotRange->{{to,tf},frange},PlotStyle->{RGBColor[0,0,0],Thickness[0.008`]},ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},DisplayFunction->Identity,Filling->{1->{2}},FillingStyle->fillcolor];p6t=Append[p6t,pinter],{k,1,nint}];p6=Show[p6t,ImageSize->{72 size,((72 size)/3)},PlotRegion->{{0,1},{0,1}},AxesLabel->{x,"error"},DisplayFunction->Identity];accumerror=Table[{to+(k-1) h,Fdiscvalues[[k,2]]-F[to+(k-1) h]},{k,1,nint}];p7=ListPlot[accumerror,PlotStyle->{PointSize[0.02`]},AxesLabel->{x,f},DisplayFunction->Identity];Show[GraphicsRow[{p2,p6}],ImageSize->{72 size,((72 size)/2)},DisplayFunction->$DisplayFunction]];

Introduction

OBJECTIVE: To visualize the process of using Riemann sums, the trapezoid rule, and Simpson's rule for approximating definite integrals, and to understand the error associated with each method.

To evaluate a definite integral, it is often necessary to use a numerical estimate of the integral in place of an exact value. This occurs in three instances: 1) when there is no simple formula for an antiderivative of the integrand; 2) when the antiderivative of the integrand is difficult to determine and/or evaluate; and, 3) when the integrand function is represented by a table of numeric values rather than by a formula. In this module, we investigate several numerical methods for estimating definite integrals.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Riemann Sums and Errors

The most direct method for estimating a definite integral is to replace it with a Riemann sum. That is, we estimate the integral as follows: (∫_a^b)f(x)dx≈(∑_(i=1)^n)f((Subscript[c, i]))h, where the interval [a, b] is divided into n subintervals, each of length h=((b-a)/n), and where (Subscript[c, i]) is any value of x taken from the (i^th) subinterval. To assess the error in the approximation, we will perform a numerical experiment wherein we use a Riemann sum to estimate an integral for which we know the exact value. (You should keep in mind, however, that numerical integration is primarily used for estimating integrals in situations where we can't determine a decimal or fraction representation for the exact value of the integral.) Our first goal is to determine what factors affect the error when we use left- or right-hand Riemann sums to estimate the exact value of an integral.

The integral we choose for our experiment is (∫_0^(π/2))cos xdx=1. First, we look at a graphical depiction of the situation. For this, we include a specially designed command that shows the areas of the rectangles that are accumulated in the Riemann sum and the errors in the estimate. The command is riemannsum[f_, x_, a_, b_, n_, rightleft_]. The arguments are the function, f, the independent variable, x, the lower bound on the integral, a, the upper bound, b, the number of rectangles, n, and a "right" or "left" sum indicator. Red areas are positive, and blue areas are negative. The error is taken as the Riemann estimate minus the exact value of the integral.

The special command riemannsum[ ] is only available in this module and is not a built-in Mathematica command. Here's how it works.

n=5;
riemannsum[Cos[x],x,0,Pi/2,n,"right"]
riemannsum[Cos[x],x,0,Pi/2,n,"left"]

The area of each small triangular-shaped region on the error graph is the local error associated with each rectangle in the Riemann sum. The sum of the areas of these triangular-shaped regions is the global or total error in the estimate of the integral. 

In general, the areas can be positive, negative, or 0. In the graphs generated by riemannsum[ ], positive areas are red, and negative areas are blue. Because the areas are signed, we will refer to them hereafter as "signed areas."

You Try It: Visualizing the Errors

Use the riemansumm[ ] command in the preceding cell (copied below) to respond to the items that follow.

n=5;
riemannsum[Cos[x],x,0,Pi/2,n,"right"]
riemannsum[Cos[x],x,0,Pi/2,n,"left"]

a)  Increase the number of rectangles, highlighted in red, in the preceding input cell, and describe qualitatively what happens to the error in the estimate of the integral each time you double the number of rectangles.

b)  What effect does the slope of the function have on the error?

c)  The function cos x decreases on the interval from 0 to (π/2), and its first derivative is negative. In this case, the right-hand Riemann sum underestimates the integral (i.e., the error is negative), and the left-hand sum overestimates the integral (i.e., the error is positive). Describe what happens to the error when we use left- and right-hand Riemann sums to estimate (∫_0^(π/2))sin xdx.

d)  While the error for each rectangle (i.e., the local error) on the left-hand sum appears to be nearly equal in magnitude to the local error for each corresponding rectangle on the right-hand sum, they aren't exactly the same because the graph of the function curves. Where does the difference in the local errors appear to be largest, and how is the difference in the local errors related to the second derivative of the function? (To see this effect more dramatically, try using two and/or three rectangles in the riemannsum[ ] command.)

e)  Based on your observations in part (d), specify two different ways in which you could use Riemann sums to improve the estimate of the integral by reducing the errors. The trick is to look for ways to estimate the integrals so that the local errors nearly cancel. For your improved methods, which feature of the integrand function would you expect to have the most significant effect on the error?

f)  Try the riemannsum[ ] command on a function that we can't integrate exactly. Two examples are (∫_1^2)(1/x)dx and (∫_0^4)(e^(-(x^2)))dx.

Part II: Analyzing the Errors

Now that we have made some qualitative observations about Riemann estimates of an integral and the associated errors, we are ready to be more quantitative in our analysis. For this we write the following riemSumLeft[f_,x_,a_,b_,n_] and riemSumRight[f_,x_,a_,b_,n_] commands to calculate left-hand and right-hand Riemann sums. The arguments are f, the integrand function, x, the independent variable, a, the lower bound, b, the upper bound, and n, the number of subintervals.

riemSumLeft[f_,x_,a_,b_,n_]:=Sum[(f/.x->(a+(i-1)*(b-a)/n))*(b-a)/n,{i,1,n}]

riemSumRight[f_,x_,a_,b_,n_]:=Sum[(f/.x->(a+i*(b-a)/n))*(b-a)/n,{i,1,n}]

Note that h=(b-a)/n and, for left-hand sums, (Subscript[c, i])=a+(i-1)h=a+(i-1)(b-a)/n, whereas for right-hand sums, (Subscript[c, i])=a+i h=a+i(b-a)/n. Inside the Sum[ ] command for the left-hand sum, the code f/.x->(a+(i-1)*(b-a)/n)) means that the value of f should be used with x replaced with (a+(i-1)*(b-a)/n)) which is (Subscript[c, i]). Similarly, x is replaced with (a+i*(b-a)/n))*(b-a)/n for the right-hand sum.

Now we use these commands, together with the fact that the exact value of the integral (∫_0^(π/2))cos xdx is 1, to build a table of values that includes the step size h, the numerical estimates of the integral, RleftI and RrightI, and the errors in the estimates, Eleft=RleftI-1 and Eright=RrightI-1. We start with two rectangles and double the number of rectangles ten times.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue];

t1=Table[{2^i,π/2/2^i,RleftI=riemSumLeft[f[x],x,a,b,2^i]//N,RrightI=riemSumRight[f[x],x,a,b,2^i]//N,RleftI-exactvalue,RrightI-exactvalue},{i,1,11}];

TableForm[t1,TableHeadings->{None,{"n","h","left est.","right est.","left error","right error"}}]

The data in the table should confirm your qualitative observations from "You Try It: Visualizing Errors" above, but now we can be more specific. We make the following quantitative observations.

1)  As the number of rectangles increases, the error decreases, and, in fact, each time we double the number of rectangles or cut h in half, we cut the error roughly in half. A numerical estimate that exhibits this characteristic is said to have an error of order h, and we designate this as O(h). 

2)  The left-hand sum overestimates the integral, whereas the right-hand some underestimates it, and, for the same number of rectangles, the errors are roughly equal in magnitude.

3)  When the number of rectangles is small, the difference between the errors for the left- and right-hand estimates is larger, showing the effect of the second derivative and the concavity of the graph on the difference in errors.

The one effect that is not evident from the numeric data is that the local error is larger when the magnitude of the first derivative, that is, the magnitude of the slope of the graph, is larger. We explore this effect further in Part III below.

You Try It: On Another Function

Perform a numerical experiment like the one in Part II on the integral (∫_0^(π/2))(sin^2)xdx, and then answer the questions that follow. To help, we include the commands to generate the table of values like the one in Part II.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue];

t1=Table[{2^i,π/2/2^i,RleftI=riemSumLeft[f[x],x,a,b,2^i]//N,RrightI=riemSumRight[f[x],x,a,b,2^i]//N,RleftI-exactvalue,RrightI-exactvalue},{i,1,11}];

TableForm[t1,TableHeadings->{None,{"n","h","left est.","right est.","left error","right error"}}]

1.  What happens to the errors each time the number of rectangles is doubled? What is the order of the error for each Riemann sum estimate of the integral?

2.  Does the left-hand Riemann sum overestimate or underestimate the integral? What about the right-hand Riemann sum?

3. What effect does the second derivative of the integrand have on the difference of the errors for the left and right Riemann sums?

Part III: The Effect of f'on the Error

To quantify the effect of the slope on the error, we investigate straight-line functions, f(x)=m x, with varying slopes.  First, let's look at a graphical depiction.

n=5;
m=2;
riemannsum[m*x,x,0,5,n,"right"]
riemannsum[m*x,x,0,5,n,"left"]

Now we set up a table of numeric values, but in this case we keep the number of rectangles fixed at 100 and change m, the slope of the straight-line graph, starting with m=0.1 and doubling it ten times. The exact values of the integrals are ((25m)/2).

Clear[m];
m[i_]=0.1*2^i;
f[x_]=m[i]  x;
a=0;
b=5;

exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]

t2=Table[{m[i],exactvalue,RleftI=riemSumLeft[f[x],x,a,b,100]//N,RrightI=riemSumRight[f[x],x,a,b,100]//N,RleftI-exactvalue,RrightI-exactvalue},{i,0,10}];

TableForm[t2,TableHeadings->{None,{"slope","exact","left est.","right est.","left error","right error"}}]

What happens to the error each time you double the slope of the line? 

Note that for nonlinear functions, the effect of the slope on the error is local in that it varies from point to point on the graph. But if we double all of the slopes, we would also double the error in using a left- or right-hand Riemann sum to estimate the integral.

Note also that the differences in the magnitudes of the errors are all 0. This supports the observation that the difference in the errors is dependent upon the second derivative or the concavity of the graph. One method for improving our numerical estimate of the integral is to add the left- and right-hand estimates and then divide the result by two. The idea that motivates this is that the errors for the two estimates will nearly cancel each other out when we add the left- and right-hand sums together. For linear functions, the second derivative is zero, and the graph has no concavity; therefore, this improvement will give the exact value of the integral.

You Try It: The Effect of f''on the Error Difference

Design a numerical experiment like the one in Part III to demonstrate the effect of the value of the second derivative on the difference between the errors in the estimate of the integral that are obtained using right- and left-hand Riemann sums. To do this, select a simple function that has a constant second derivative, is easy to integrate, and allows you to change the value of the second derivative. Generate a table of numerical values that shows the values of the second derivative and the difference in the errors for the left and right Riemann sums. Take 1 as the starting value of the second derivative, and double its value ten times. Try functions that are concave down as well as ones that are concave up, and summarize the results of your investigation. (We have put the solution for this problem in the next input cell with the code hidden. Try to do it on your own, and then look at our solution. To reveal the hidden code, select the small cell bracket to the right, pull down the Cell menu, select Cell Properties, and click on Cell Open.)

Clear[ac];
ac[i_]=2^i;
f[x_]=(1/2)ac[i] (x^2);
a=0;
b=5;
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue];

t2=Table[{ac[i],exactvalue//N,RleftI=riemSumLeft[f[x],x,a,b,100]//N,RrightI=riemSumRight[f[x],x,a,b,100]//N,RleftI-exactvalue,RrightI-exactvalue,(RleftI-exactvalue)-(RrightI-exactvalue)},{i,0,10}];

TableForm[t2,TableHeadings->{None,{"f''(x)","exact","left estimate","right estimate","left error","right error","error difference"}}]

Part IV: Trapezoids and Midpoint Riemann Sums

In this part, we investigate two ways to improve the efficiency of numerical estimates of definite integrals.

Trapezoids

We achieve the first improvement by adding together the left and right Riemann sum estimates of the integral and dividing by two. The motivation for this approach is that the errors in the left and right sums should nearly cancel each other out. In the next cell, we form the new function trapezoid[ ] that does this.

trapezoid[f_,x_,a_,b_,n_]:=(riemSumLeft[f,x,a,b,n]+riemSumRight[f,x,a,b,n])/2

(The trapezoid method gets its name from the signed area of a trapezoid in each subinterval that results from adding the signed areas of the left and right rectangles in each subinterval and dividing by 2. That is, (f((Subscript[c, left]))h+f((Subscript[c, right]))h)/2=((f((Subscript[c, left]))+f((Subscript[c, right])))/2)h, the area of a trapezoid, where f((Subscript[c, left])) is the length of one parallel side, f((Subscript[c, right])) is the length of the other parallel side, and h is the distance between the two sides.)

Now we use trapezoid[ ] to estimate (∫_0^(π/2))cos xdx with ten subintervals and compare the result with the left- and right-hand Riemann sums. First, we calculate the exact value of the integral.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]

Now we calculate the estimates of the integral.

n=10;
trapezoid[f[x],x,a,b,n]//N
riemSumLeft[f[x],x,a,b,n]//N
riemSumRight[f[x],x,a,b,n]//N

And we compute the errors.

(trapezoid[f[x],x,a,b,n]//N)-exactvalue
(riemSumLeft[f[x],x,a,b,n]//N)-exactvalue
(riemSumRight[f[x],x,a,b,n]//N)-exactvalue

For the same number of subintervals, we get a much more precise estimate of the integral from the trapezoid method.

Midpoint Riemann Sums

The second improvement we consider is to use the midpoint of each subinterval in a Riemann sum. The effect of doing this is to nearly cancel the error on either side of each midpoint in each subinterval. If the function we integrate is increasing over a subinterval, then the midpoint rectangle will overestimate the actual area on the left side of the midpoint and will underestimate, by nearly the same amount, the actual area on the right side of the midpoint. Again we expect that these errors will nearly cancel each other out. Here is the midpoint[ ] command.

riemSumMid[f_,x_,a_,b_,n_]:=Sum[(f/.x->(a+(i-1/2)*(b-a)/n))*(b-a)/n,{i,1,n}]

Let's try it on our old friend and compare the result with left and right Riemann sums and with the trapezoid estimate. First, we calculate the exact value of the integral.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]

Now we calculate the estimates of the integral.

n=10;
riemSumMid[f[x],x,a,b,n]//N
trapezoid[f[x],x,a,b,n]//N
riemSumLeft[f[x],x,a,b,n]//N
riemSumRight[f[x],x,a,b,n]//N

Again, we compute the errors.

(riemSumMid[f[x],x,a,b,n]//N)-exactvalue
(trapezoid[f[x],x,a,b,n]//N)-exactvalue
(riemSumLeft[f[x],x,a,b,n]//N)-exactvalue
(riemSumRight[f[x],x,a,b,n]//N)-exactvalue

For the same number of intervals, we get a much more precise estimate of the integral using the midpoint rule; it is even more precise than the estimate given by the trapezoid rule. What relationship does the error for the midpoint estimate seem to have with that for the trapezoid estimate?

You Try It: On Some Other Functions

Compare the Riemann left, Riemann right, trapezoid, and midpoint estimates of some integrals that you pick. How does the error for the midpoint rule compare with the error for the trapezoid rule with the same number of subdivisions? Also, try increasing the number of subintervals, highlighted in blue. What happens to the errors when you double the number of subintervals? 

Here are some integrals that you might want to try: (∫_0^(π/2))sin xdx, (∫_(-1)^1)(e^x)dx, (∫_0^1)(Sqrt[1-(x^2)])dx. To help you out, we provide some commands in the two cells that follow.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
n=10;

Print["The exact value of the integral is ",exactvalue]
Print["The midpoint estimate is ",riemSumMid[f[x],x,a,b,n]//N]
Print["The trapezoid estimate is ",trapezoid[f[x],x,a,b,n]//N]
Print["The left Riemann estimate is ",riemSumLeft[f[x],x,a,b,n]//N]
Print["The right Riemann estimate is ",riemSumRight[f[x],x,a,b,n]//N]

Print["The midpoint error is ",(riemSumMid[f[x],x,a,b,n]-exactvalue)//N]
Print["The trapezoid error is ",(trapezoid[f[x],x,a,b,n]-exactvalue)//N]
Print["The left Riemann error is ",(riemSumLeft[f[x],x,a,b,n]-exactvalue)//N]
Print["The right Riemann error is ",(riemSumRight[f[x],x,a,b,n]-exactvalue)//N]

Part V: Experiment with Trapezoids and Midpoints

In this part, we do some numerical experiments to determine what factors affect the errors when we use trapezoids and midpoint Riemann sums to estimate a definite integral.

The Effect of the Step Size on the Error

First, we would like to determine the effect of the step size h on the error. We build a table similar to the one in Part II.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]

t1=Table[{2^i,π/2/2^i,trap=trapezoid[f[x],x,a,b,2^i]//N,mid=riemSumMid[f[x],x,a,b,2^i]//N,trap-exactvalue,mid-exactvalue},{i,1,11}];

TableForm[t1,TableHeadings->{None,{"n","h","trap est.","mid est.","trap error","mid error"}}]

The Effect of f'' on the Error

Based on the observations made in Parts II, III, and IV, we conjecture that the error for these methods depends on the magnitude of the second derivative of the integrand function. We test our conjecture by comparing the error in estimating the area under the quadratic functions of the form f(x)=((a (x^2))/2) between x=0 and x=5 for various values of a. The exact values of the integrals are (∫_0^5)((a (x^2))/2)dx=((125a)/6), and the second derivative of f(x) is equal to a. We set up a table in which a starts at 0.1 and doubles ten times. The number of subintervals is kept constant at 100.

Clear[ac];
ac[i_]=0.1*2^i;
f[x_]=((ac[i] (x^2))/2);
a=0;
b=5;
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]

t2=Table[{ac[i],exactvalue,trap=trapezoid[f[x],x,a,b,100]//N,mid=riemSumMid[f[x],x,a,b,100]//N,trap-exactvalue,mid-exactvalue},{i,0,10}];

TableForm[t2,TableHeadings->{None,{"f''(x)","exact","trap est.","mid est.","trap error","mid error"}}]

In "You Try It" below, we ask you to summarize the results of these experiments.

You Try It: Summarize the Experiment's Results

Summarize the results of the numerical experiments in Part V. In your summary, address the following items.

a)  How is the error for each method related to the number of subintervals n and the interval width h? What is the order of the error for each method? (Note: If cutting the interval width in half reduces the error by a factor of ((((1/2)))^n), then the error is of order (h^n), and we designate this by O((h^n)).)

b)  How is the error for the trapezoid method related to the error for the midpoint Riemann sum?

c)  How is the error related to the second derivative of the integrand function? (To answer this completely, repeat the second experiment in Part V for negative values of a. 

d)  Are your conclusions consistent with the trapezoid error formula that is found in the text?

e)  From the pattern of errors for the trapezoid method and the midpoint Riemann sum method, specify how you might combine these two estimates of the integral so that the errors will nearly cancel each other. Compare your idea for combining the estimates with those of other students.

Part VI: Simpson's Method

Noting that the magnitude of the error for the trapezoid method is approximately two times the magnitude of the error for a midpoint Riemann sum with the same number of subintervals, we can make yet another improvement in the efficiency of our numerical estimations of the integrals. If we add two midpoint estimates and one trapezoid estimate and divide the result by three, the errors for the two methods should nearly cancel. The result of this combination is called Simpson's method.

To illustrate this, consider the following estimates of the integral (∫_0^(π/2))cos xdx. Using 10 subintervals, we compare the midpoint and trapezoid estimates and their errors, and then we combine them as indicated above and calculate the error for the new estimate.

f[x_]=Cos[x];
a=0;
b=(π/2);
n=10;
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
trap=trapezoid[f[x],x,a,b,n]//N;
mid=riemSumMid[f[x],x,a,b,n]//N;

simps=(2*mid+trap)/3;

Print["The exact value of the integral is ",exactvalue]
Print["The trapezoid estimate is ",trap]
Print["The trapezoid error is ",(trap-exactvalue)//N]
Print["The midpoint estimate is ",mid]
Print["The midpoint error is ",(mid-exactvalue)//N]
Print["The Simpson estimate of the integral is ",simps]
Print["The error for Simpson's estimate is ",simps-exactvalue]

Since Simpson's method requires three function evaluations in each subinterval, one at each end point for the trapezoid plus one at the mid point, we count the number of subintervals differently. If we take the dividing points for the subintervals to be all of the places where the integrand function is evaluated, then there are actually twice as many subintervals for Simpson's method than there would be for the trapezoid method or the midpoint Riemann sum. One consequence of this is that we need an even number of subintervals for Simpson's method to work. 

In the next cell we write a command that gives the Simpson estimate of a definite integral.

simpson[f_,x_,a_,b_,n_]:=(2*riemSumMid[f,x,a,b,n/2]+trapezoid[f,x,a,b,n/2])/3

Let's use it to estimate (∫_0^(π/2))cos xdx=1 and then calculate the errors, comparing the results with those obtained using trapezoids and midpoint Riemann sums.

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue]
Print["The Simpson estimate of the integral is ",simpson[Cos[x],x,a,b,6]//N]
Print["The trapezoid estimate of the integral is ",trapezoid[Cos[x],x,a,b,6]//N]
Print["The midpoint estimate of the integral is ",riemSumMid[Cos[x],x,a,b,6]//N]

Now we compare the errors.

Print["The error in the Simpson estimate is ",(simpson[Cos[x],x,a,b,6]//N)-exactvalue]
Print["The error in the trapezoid estimate is ",(trapezoid[Cos[x],x,a,b,6]//N)-exactvalue]
Print["The error in the midpoint estimate is ",(riemSumMid[Cos[x],x,a,b,6]//N)-exactvalue]

Using only six subintervals, Simpson's method gives very good results.

While left and right Riemann sums integrate constant functions exactly, and the trapezoids and midpoint Riemann sums integrate linear functions exactly, we might conjecture that Simpson's rule integrates quadratic functions exactly. This conjecture is in fact correct, but it is even better than that. Simpson's method integrates cubic functions exactly. We show this in the following cell by applying Simpson's method to a general cubic function with only two subintervals and then comparing the result with the exact value.

Clear[fcubic,a,b];
fcubic[x_]=c (x^3)+d (x^2)+e x +f;
simpsonresult=simpson[fcubic[x],x,a,b,2]//Expand

integralresult=Integrate[fcubic[x],{x,a,b}]//Expand

simpsonresult==integralresult

You Try It: Experiment with Simpson

Perform some numerical experiments like those in the preceding parts to show that the order of the error for Simpson's method is fourth. Also show that the local error is proportional to the fourth derivative of the integrand function. (We have put the solutions for this problem in the next two input cells with the code hidden. Try to do it on your own, and then look at our solutions. To reveal the hidden code, select the small cell bracket to the right of each cell, pull down the Cell menu, select Cell Properties, and click on Cell Open.)

f[x_]=Cos[x];
a=0;
b=(π/2);
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue];

t1=Table[{2^i,π/2/2^i,simps=N[simpson[f[x],x,a,b,2^i],20],simps-exactvalue},{i,1,11}];

TableForm[t1,TableHeadings->{None,{"n","h","Simpson estimate","Simpson error"}}]

Clear[ac];
ac[i_]=0.1*2^i;
f[x_]=((ac[i] (x^4))/(4!));
a=0;
b=5;
exactvalue=(∫_a^b)f[x]\[DifferentialD]x;
Print["The exact value of the integral is ",exactvalue];

t2=Table[{ac[i],SetPrecision[exactvalue,12],simps=SetPrecision[simpson[f[x],x,a,b,100],12],simps-exactvalue},{i,0,10}];

TableForm[t2,TableHeadings->{None,{"(f^((4)))(x)","exact","Simpson estimate","Simpson error"}}]