(C) Copyright 1994 D.C.M. Burbulla and  C.T.J. Dodson
       Notebook to support the book
     Self-Tutor for Computer  Calculus

Using Mathematica 2

by
     D. C. M. Burbulla and C. T. J. Dodson, University of Toronto
 Published by:    
    Prentice-Hall Canada  (English)   ISBN 0130152803     and 
     Prentice-Hall Japan  (Japanese) ISBN4-8101-8558-3 C3055 P2880E

Numerical Integration

As we have seen in the previous chapter, Riemann sums are not a practical
way to calculate definite integrals: the Fundamental Theorem of Calculus 
provides a much more direct method of evaluating a definite integral, if an 
antiderivative of the integrand is known. However, it is amathematical fact of 
life that not all elementary functions have elementaryantiderivatives. Thus there 
is still a need for numerical methods to approximate some definite integrals.
In this chapter we introduce some functions to augment those  built into Mathematica.

Chap8.m

riemann/: riemann[f_, {a_, b_}, n_, choice_:1] :=
    (b - a)*Sum[f[(1 - choice)*(a + (i - 1)*(b - a)/n) + 
    choice*(a + i*(b - a)/n)], {i, 1, n}]/n // N

riemann::usage = "riemann[f, {a, b}, n, choice] calculates the numerical 
value of the Riemann sum obtained by dividing the interval {a, b} 
into n equally spaced subintervals and evaluating f at the chosen point of each 
subinterval. The chosen point is determined by a number between 0 and 1, ranging 
from the left end point to the right end point of each subinterval. The default 
choice is 1, ie. choosing the right end point of each subinterval."


trapezoid/: trapezoid[f_, {a_, b_}, n_] := 
 ((b - a)*(f[a] + 2*Sum[f[a + (i*(b - a))/n], {i, 1, n - 1}] +
                    f[b]))/(2*n) // N
trapezoid::usage = 
"trapezoid[f, {a, b}, n] calculates the trapezoid rule
approximation of the definite integral of f on the interval {a, b} by
dividing the interval {a, b} into n equally spaced subintervals."

  
simpson/: simpson[f_, {a_, b_}, n_] := 
   ((b - a)*(f[a] + 4*Sum[f[a + (i*(b - a))/n], {i, 1, n - 1, 2}] + 
 2*Sum[f[a + (j*(b - a))/n], {j, 2, n - 2, 2}] + f[b]))/(3*n) // N

simpson::usage = 
"simpson[f, {a, b}, n] calculates the Simpson's rule
approximation of the definite integral of f on the interval {a, b} by
dividing the interval {a, b} into n equally spaced subintervals, 
where n is even."


er/: er[f_, {a_, b_}, n_, choice_:1] :=
   Abs[NIntegrate[f[x], {x, a, b}] - riemann[f, {a, b}, n, choice]]

er::usage = "er[f, {a, b}, n, choice] computes the absolute value of the
difference between the value of the definite integral of f on the 
interval {a, b} and the value of the Riemann sum approximation 
riemann[f, {a, b}, n, choice]."


et/: et[f_, {a_, b_}, n_] := 
   Abs[NIntegrate[f[x], {x, a, b}] - trapezoid[f, {a, b}, n]]
 
et::usage = "et[f, {a, b}, n] computes the absolute value of the
difference between the value of the definite integral of f on the 
interval {a, b} and the value of the trapezoid rule approximation 
trapdezoid[f, {a,b}, n]." 


es/: es[f_, {a_, b_}, n_] :=
   Abs[NIntegrate[f[x], {x, a, b}] - simpson[f, {a, b}, n]]
 
es::usage = "es[f, {a, b}, n] computes the absolute value of the
difference between the value of the definite integral of f on the 
interval {a, b} and the value of the Simpson's rule approximation 
simpson[f, {a,b}, n]." 


s/: s[f_, {a_, b_}, x_, choice_:1] := 
	If[a < x < b, f[(1 - choice)*a + choice*b], 0]
 
 
q/: q[f_, {a_, b_}, x_] := 
   If[a < x < b, f[a] + (2*(f[(a + b)/2] - f[a])*(x - a))/(b - a) + 
(2*((f[a] + f[b]) - 2*f[(a + b)/2])*(x - a)*(x - (a + b)/2))/
                    (b - a)^2, 0]

 
l/: l[f_, {a_, b_}, x_] := 
   If[a < x < b, ((f[b] - f[a])*(x - a))/(b - a) + f[a], 0]

 
viewrapprox/: viewrapprox[f_, {a_, b_}, n_, choice_:1] := 
 (one = Block[{x}, Plot[f[x], {x, a - 0.1, b + 0.1},
  DisplayFunction -> Identity]];
     two = Block[{i, x}, 
Plot[Release[Table[s[f, {a + ((i - 1)*(b - a))/n, 
a + (i*(b - a))/n}, x, choice], 
      {i, n}]], {x, a - 0.1, b + 0.1}, PlotRange -> All, PlotPoints -> 50,
            DisplayFunction -> Identity]];
      Block[{x}, Show[two, one, PlotLabel -> 
     {f[x], n*"intervals", N[(b - a)/n]*"= norm", 
If[Chop[er[f, {a, b}, n, choice]] == 0, "no error", 
er[f, {a, b}, n, choice]*"= error"]},
 PlotRange -> All, 
 DisplayFunction -> $DisplayFunction ]])
 
viewrapprox::usage = 
"viewrapprox[f, {a, b}, n, choice] shows the Riemann 
sum approximation to the definite integral of f on the interval 
{a, b} for n equally spaced subintervals of the given interval 
{a, b}. The default choice in each subinterval is the right end point. 
The PlotLabel includes the function f, the number of 
subintervals, the norm of the partition, and the absolute error 
of the Riemann approximation."

 
viewtapprox/: viewtapprox[f_, {a_, b_}, n_] := 
(one = Block[{x}, Plot[f[x], 
{x, a, b}, DisplayFunction -> Identity]];
      two = Block[{x, i}, Plot[Release[Table[
     l[f, {a + ((i - 1)*(b - a))/n, a + (i*(b - a))/n}, x],
      {i, n}]], {x, a, b}, PlotRange -> All, PlotPoints -> 50,
            DisplayFunction -> Identity]];
      Block[{x}, Show[two, one, PlotLabel ->
     {f[x], n*"intervals", N[(b - a)/n]*"= norm",
If[Chop[et[f, {a, b}, n]] == 0, "no error", 
et[f, {a, b}, n]*"= error"]}, PlotRange -> All, 
DisplayFunction -> $DisplayFunction ]])

viewtapprox::usage = 
"viewtapprox[f, {a, b}, n] shows the trapezoid 
rule approximation to the definite integral of f on 
the interval {a, b} 
for n equally spaced subintervals of the given interval {a, b}. 
The PlotLabel includes the function f, the number of subintervals, 
the norm of the partition, and the absolute error of the 
trapezoid rule approximation."


viewsapprox/: viewsapprox[f_, {a_, b_}, n_] := 
 (one = Block[{x}, Plot[f[x], {x, a, b}, 
DisplayFunction -> Identity]]; two = Block[{x, i},
      Plot[Release[Table[
   q[f, {a + ((i - 1)*(b - a))/n, a + ((i + 1)*(b - a))/n}, x],
      {i, 1, n - 1, 2}]], {x, a, b}, 
      PlotRange -> All, PlotPoints -> 50,
            DisplayFunction -> Identity]];
      Block[{x}, Show[two, one, PlotLabel ->
     {f[x], n*"intervals", N[(b - a)/n]*"= norm",
    If[Chop[es[f, {a, b}, n]] == 0, "no error",  
es[f, {a, b}, n]*"= error"]},
PlotRange -> All, DisplayFunction -> $DisplayFunction]])

viewsapprox::usage = "viewsapprox[f, {a, b}, n] 
shows the Simpson's rule approximation to the definite integral of f 
on the interval {a, b} for n equally spaced subintervals of the given 
interval {a, b}, where n must be even. 
The PlotLabel includes the function f, the number of subintervals,
the norm of the partition, and the absolute error of the 
Simpson's rule approximation."

?riemann

The Trapezoid Rule and Simpson's Rule

?trapezoid

?simpson

Error expressions are included as well:

?er

?et

?es

Examples

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

trapezoid[f, {0, 1}, 4]

You can view graphs of the function f and your trapezoid rule
approximation by using viewtapprox. Similar functions are defined to view 
Riemann sum and Simpson's rule approximations

viewtapprox[f, {0, 1}, 4];

trapezoid[f, {0, 1}, 10]

simpson[f, {0, 1}, 10]

viewsapprox[f, {0, 1}, 10]

We can see how good Simpson's Rule is by looking at the reciprocal of 
the error:

Table[{n, 1/es[f, {0, 1}, n]}, {n, 2, 10, 2}]

ListPlot[%];

Try similar error plots using et and er.