Basic Approximation
of
Definite Integrals

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

The Functions leftRect[.],  midRect[.], trapezoids[.]:

Off[General::spell1];
Needs["Graphics`FilledPlot`"];
Needs["Graphics`Colors`"];

leftRect[n_,f_,a_,b_] :=
Module[ {X,dX,rect,rects,partition,fctplot,leftSum},
fctplot = Plot[ f[x], {x, a, b},
		DisplayFunction -> Identity ];
	rect[x_,y_] := Line[{{x,0},{y,0},{y,f[x]},{x,f[x]},
		{x,0}}];
	dX = (b - a) / n ;
	X = Table[ a + i (b - a) / n ,{i, 0, n}];
	partition := Graphics[ {GrayLevel[0.7],
			Table[ Rectangle[ {X[[i]], 0}, {X[[i+1]], 
				f[X[[i]]]} ], {i, 1, n } ] } ] ;
	rects := Graphics[Table[rect[X[[j]],X[[j+1]]],
			{j,1,n}]];
	leftSum = N[ Sum[ f[X[[i]]] dX, {i, 1, n } ], 10];
	label = leftSum "= area of rectangles "  ;
	Show[partition, rects, fctplot,PlotLabel -> label,
	DisplayFunction->$DisplayFunction]; 
	Return[leftSum]
	] ;

midRect[n_,f_,a_,b_] :=
Module[ {X,dX,mid,mids,partition,fctplot,midSum},
fctplot = Plot[ f[x], {x, a, b},
		DisplayFunction -> Identity ];
		mid[x_,d_] := Line[{{x,0}, {x+d,0},
		{x+d,f[x + d/2]}, {x,f[x + d/2]}, {x,0}}] ;
		dX = (b - a) / n ;
		X = Table[ a + i (b - a) / n ,{i, 0, n}];

		partition := Graphics[ {GrayLevel[0.7],
				Table[ Rectangle[ {X[[i]], 0}, {X[[i+1]], f[X[[i]]+dX/2]} ],
				 {i, 1, n} ] } ] ;
		mids := Graphics[Table[mid[X[[j]],dX],
				{j,1, n } ] ] ;
		midSum = N[ Sum[ f[X[[i]]+dX/2] dX, {i,1, n} ], 10];
		label = midSum "= area of rectangles" ;
		Show[ partition, mids, fctplot, PlotLabel -> label,
		DisplayFunction->$DisplayFunction ]; 
		Return[midSum]
	];

trapezoids[n_,f_,a_,b_]:=
Module[ {X,dX,trap,trapz,partition,fctplot,trapSum},
fctplot = Plot[ f[x], {x, a, b},
		DisplayFunction -> Identity ];
		dX = (b - a) / n ;
		X = Table[a + i (b-a)/n,{i,0,n}];
		trap[x_,d_]:=Line[{{x,0},{x+d,0},{x+d,f[x+d]},
			{x,f[x]},{x,0}}];

		partition := Graphics[ {GrayLevel[0.7],
			Table[Polygon[{{X[[i]], 0},{X[[i]],f[X[[i]]]},
			{X[[i+1]], f[X[[i+1]]]},{X[[i+1]],0}}],
				 {i, 1, n} ] } ] ;
		trapz := Graphics[Table[trap[X[[j]],dX],
			     {j, 1, n} ]];		 
		
		trapSum = N[dX (f[a] + 2 Sum[f[X[[i]]],{i,2, n }] 
			+ f[b])/2, 10];
		label = trapSum"= area of trapezoids"  ;
		Show[ partition, trapz, fctplot, PlotLabel -> label, DisplayFunction->$DisplayFunction ]; 
		Return[trapSum]
	];

Overview: Approximating Area Under a Curve

A Riemann sum approximates the area under a curve by dividing the region into rectangles.  The sum of the areas of the rectangles is a better and better approximation of the actual area under the curve if the rectangles are taken smaller and smaller.

The function

leftRect[n]

below divides the area under a function curve into n rectangles and then sums the areas of these rectangles to approximate the area under the curve.  By activating the function cell and typing leftrect[n], Mathematica will partition the area under the function f[x] = 1/(1 + x^2)  into  n  rectangles and sum over this partition.  Mathematica also draws a graph of your partition and gives you the true value so you can see how closely the sum approximates the actual area.  You may change the function and limits by changing the definition of f[x] or a and b below.

Approximation of area does not require that we take the approximating rectangles at the left side; we can also use rectangles that touch the curve on the right or in the middle.  Midpoint rectangles tend to be the most accurate.

We do not even need to use rectangles (with a horizontal top).  A trapezoid that matches at both the left and right is another good approximation, comparable to midpoint rectangles.  Parabolic regions could also be used as in "Simpson's Approximation."  The parabolic approximation is the same as the weighted average of (1/3) trapezoids + (2/3) midpoint rectangles.  However, Mathematica has its own NIntegrate[.] function that we can use for accurate numerical approximations to integrals, so the main point of this NoteBook is that the approximations converge, not how fast they converge.

References to the Text

This NoteBook goes with Chapter 12 of the text.

Step-by-step

The function and limits

We want to illustrate the sum approximation for the area interpretation of the definite integral.  We will use the function

f[x] = (4/(1 + (x^2)))

because this has an interesting answer which is not so obvious symbolically.  The first computation below here shows the area we want to approximate:

\(Clear[f, x]; 
f[x_] := 4\/\(1 + x\^2\); 
a = 0; 
b = 1; 

    FilledPlot[f[x], {x, a, b}, PlotStyle -> {{Thickness[0.02], Red}}, 
      Fills -> {Blue}]; \)

The next computation shows how to use the leftRect[.], midRect[.], and trapezoid[.] functions to both draw the approximation and calculate its value.  The input to each of these functions is the same,

leftRect[ number of subdivisions, function f[x], left endpoint  a, right endpoint b ]

<Enter> the next sample:

n = 4; leftSum = leftRect[n, f, a, b]; midSum = midRect[n, f, a, b]; trapSum = trapezoids[n, f, a, b]; Print["With ", n, " subdivisions:"]; Print["Left rectangles total =", leftSum]; Print["Midpoint rectangles total =", midSum]; Print["Trapezoids total =", trapSum]; Print["Accurate answer =", NIntegrate[f[x], {x, a, b}]];

Practice Exercise

Draw and calculate the left, mid and trapezoid approximations under the graph of

f[x] = (4/(1 + (x^2)))  from  a = 0  to  b = 1  with 10 subdivisions:

Left Rectangles

<Enter> and run the Mathematica animation in the graphics cell below. 
(Press <Command> & <Y> for animation.)

Do[leftSum = leftRect[(2^n), f, a, b], {n, 1, 6}]; Print["Actual Area (Definite Integral) = ", 
      N[(∫_a^b) f[x] \[DifferentialD]x, 10]];

MidPoint Rectangles

<Enter> and run the Mathematica animation in the graphics cell below.

Do[midSum = midRect[(2^n), f, a, b], {n, 1, 6}]; Print["Actual Area (Definite Integral) = ", 
      N[(∫_a^b) f[x] \[DifferentialD]x, 10]];

Trapezoids

<Enter> and run the Mathematica animation in the graphics cell below.

Do[trapSum = trapezoids[(2^n), f, a, b], {n, 1, 6}]; Print["Actual Area (Definite Integral) = ", 
      N[(∫_a^b) f[x] \[DifferentialD]x, 10]];

Numerical summary of the loops to (2^6):

Print["With ", (2^6), " subdivisions,"]; Print["Left Rectangles = ", leftSum]; Print["Mid Rectangles = ", midSum]; Print["Area of Trapezoids = ", trapSum]; Print["Simpson's rule = (2/3)Midpts + (1/3)Traps =", 
      ((2 midSum)/3) + (trapSum/3)]; Print["Area to more than machine precision = ", 
      N[(∫_a^b) f[x] \[DifferentialD]x, 10]]; Print["Exact symbolic area = ", 
      (∫_a^b) f[x] \[DifferentialD]x];

Practice Exercise

Use leftRect[.], midRect[.] and trapezoids[.] with lots of subdivisions to approximate the area under  y = (x^2), for  x  from  0  to  1.  You probably already know that the integral is exactly computable.  What is the area? (Based on observation or extra knowledge.)

Exercises

Exercise 12.3.1: Slow convergence to a familiar answer

Use leftRect[.], midRect[.] and trapezoids[.] with lots of subdivisions to approximate the area of the unit circle. (HINT: Define the function  f[x] = (Sqrt[1 - (x^2)]) and make a = -1, b = +1.  Multiply your answer or  f[x]  by  2 to get  Pi.)

Exercise 12.8.1: Numerical Approximation

Run the comparison computations belos.  Type answers to the questions following the computations and save the NoteBook.

We will make some numerical comparisons on the accuracy of various sum approximations to the integral of  f[x] dx  when

f[x] = (4/(1 + (x^2)))

To begin, enter the function definiton in the next cell.

f[x_] := (4/(1 + (x^2)))

Straight Summation

This computation uses the simple summation command with no refinements. (Be sure you entered the definition of  f[x].)

a = 0b = 1dx =  .1(∑_({{x = a}, {Δ x = dx}})^(b - dx))f[x] dx

dx =  .0001NSum[f[x] dx, {x, a, b - dx, dx}]

Questions:

1)  How many terms are added in the two sums above?

Your answer 1:

2)  Compute the sum with 1,000 and 5,000 terms. (Using the same  f[x], a, and b.)

3)  How accurate is the answer compared with dx?

Your answer 3:

Midpoint Summation

midsum[f_, a_, b_, dx_]:=(∑_({{x = a + (dx/2)}, {Δ x = dx}})^(b - (dx/2)))f[x] dx

a = 0b = 1dx =  .01midsum[f, a, b, dx]

dx =  .001midsum[f, a, b, dx]

Questions:

Q.1.  How many terms are added in the two sums above?

Your answer 1:

Q.2.  How accurate is the answer compared with dx?

Your answer 2:

Trapezoidal Summation

trapsum[f_, a_, b_, dx_]:=(1/2) (f[a]+2 (∑_({{x = a + dx}, {Δ x = dx}})^(b - dx))f[x]+f[b]) dx

a = 0b = 1dx =  .01trapsum[f, a, b, dx]

dx =  .001trapsum[f, a, b, dx]

Questions:

Q.1.  How many terms are added in the two sums above?

Your answer 1:

Q.2.  How accurate is the answer compared with dx?

Your answer 2:

Simpson's Summation

simpson[f_, a_, b_, dx_] := 
      (1/3) trapsum[f, a, b, dx] + (2/3) midsum[f, a, b, dx]

a = 0b = 1dx =  .01simpson[f, a, b, dx]

dx =  .1simpson[f, a, b, dx]

Questions:

Q.1.  How many terms are added in the two sums above?

Your answer 1:

Q.2.  How accurate is the answer compared with dx?

Your answer 2:

Mathematica Numerical Integration

a = 0b = 1NIntegrate[f[x], {x, a, b}]

Questions:

Q.1.  How fast did Mathematica compute the approximate integral above?

Your answer 1:

Q.2.  How accurate is the answer?

Your answer 2:

Q.3.  Did you need to tell Mathematica dx??

Your answer 3:

Mathematica Symbolic Definite Integration

a = 0b = 1(∫_a^b) f[x] \[DifferentialD]x

N[π]

Mathematica Symbolic Indefinite Integration
                   (with & without limits)

Clear[a, b]a := 0b := 1f[x_] := (1/(1 + (x^2)))(∫_a^b) f[x] \[DifferentialD]x

Integrate[f[x],x]

The Big Picture

Question 1: 
How does each approximation compare to the final limit on a scale of the increment?

For example, the third left rectangle approximation to the area under  4/(1 + (x^2))  from  0  to  1  has  8  steps, each of size  1/8.  The approximating sum is  0.815997  and the actual area is  0.785398, for an actual error of  0.815997 - 0.785398 = 0.030599.  This error is too high because the curve is decreasing, so that left rectangles extend beyond the region.

The increment of 1/8  is 0.125, so the error is  0.030599/0.125 = 0.244792, or about  1/4 of the size of the x-step.

The sixth step of the left rectangle approximation has  2^6 = 64 steps of size 1/64.  The actual error is   0.7892942409 - 0.7853981634 = 0.00389608  and the error is the fraction    0.249349  of the step size, or about  1/4 of the step size of  1/64.

It looks like the error is about 1/4 the size of the step.  In any case, the order of the size of the step.  What is the scale of the errors for the other methods.  Δ x, (((Δ x))^2)?

Question  2:
Why is the midpoint rectangle a better approximation to the area than the left rectangle?

Question  3:
How do the errors compare for the midpoint rule and the trapezoid rule?

Beyond the Basic Theory: Riemann Sums (Optional)

Random Riemann Sums (Optional)

A general Riemann sum (needed for discontinuous functions) can have variable size partitions and random points of evaluation.  As the number of rectangles approaches infinity, the Riemann sum will still approach the area under the curve provided that the size of the biggest increment in  x  tends to zero.  When you activate the following function cell and type riemannsum[n], Mathematica will divide the area under the curve into rectangles of random widths and sum their areas.

Riemann Function - <Enter> this, but don't read the code.

f[x_] := 1/(1 + x^2);
a = 0 ;
b = 1 ;
riemannsum[n_] :=
Module[ {dx},
		RandRect[f_,x1_,x2_,y_]:= Line[{{x1,0},{x2,0},
			{x2,f[y]},{x1,f[y]},{x1,0}}];
		maxstep := 2 (b - a) / n ;
		x[0] = a ;
		x[n] = b ;
		Do[ 
			upperlimit = Min[ x[i-1] + maxstep, b ] ;
			x[i] = Random[ Real, {x[i-1], upperlimit} ] ;
			xp[i] = Random[ Real, {x[i-1], x[i]} ] ;
			dx[i] = x[i] - x[i-1],
			{i,1,n-1}
		  	] ;
		xp[n] = Random[ Real, {x[n-1], x[n]} ] ;
		dx[n] = x[n] - x[n-1] ;
		functionplot := Plot[ f[x], {x, a, b}, DisplayFunction -> Identity ] ;
		partition := Graphics[ { GrayLevel[0.4],
			Table[ Rectangle[ {x[i-1], 0},
				{x[i], f[ xp[i] ]} ], {i, 1, n} ] } ] ;
		randrectangles:=Graphics[Table[RandRect[f,x[i-1],
			x[i],xp[i]], {i,1,n}]  ];
		Show[ partition, randrectangles, functionplot, DisplayFunction -> $DisplayFunction ] ;
		Print[ "Riemann Sum = ",
			N[ Sum[ f[xp[i]] dx[i], {i, 1, n} ], 10 ] ] ;
		Print[ "Actual Area (Definite Integral) = ",
			N[ Integrate[ f[x], {x, a, b} ], 10 ] ] ;
	] ;

A Sample Riemann Sum

riemannsum[20]