Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 7.	Numerical Integration
	Algorithm 7.1	Composite Trapezoidal Rule	
	Algorithm 7.2	Composite Simpson Rule	
	Algorithm 7.3	Recursive Trapezoidal Rule	
	Algorithm 7.4	Romberg Integration	
	Algorithm 7.5	Adaptive Quadrature Using Simpson's Rule	
	Algorithm 7.6	Gauss-Legendre Quadrature

Algorithm 7.1  (Composite Trapezoidal Rule). To approximate the integral
aòbÊf(x)dx  ~  h [f(a) + f(b)]/2  +  h Sk=1k=m-1Êf(xk)Ê
by sampling f(x) at the m+1 equally spaced points  xk = a + hk  
for  k = 0, 1, 2,...,m.   Notice that  x0 = a  and  xm = b.
Section 7.2 Composite Trapezoidal & Simpson's Rule   Page 365

Off[General::spell1];
Clear[TrapezoidalRule,f0,a,b,m];
TrapezoidalRule[f0_,a_,b_,m_] :=
Module[{f,k},
  Set@@{f[x_],f0};
  h = (b - a)/m//N;
  x0 = a//N;
  s = 0;
  For[k=1,k<=m-1,k++,
    Module[{},
      xk = x0 + h k;
      s = s + f[xk]//N; ]];
  s = h ( f[a] + f[b] )/2 + h s//N;
  Trule = s;
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to the integral over"];
  Print["[",a,",",b,"] computed by the trapezoidal rule"];
  Print["using ",m," subintervals is   ",N[Trule,14]]; ];
On[General::spell1];

Chapter 7, Example 7.5, Page 358.

Integrate  f(x) = 2 + sin( 2 x^1/2 )  over the interval [1,6].

Clear[f,x];
f[x_] = 2 + Sin[2Sqrt[x]]

line1 = Graphics[{Line[{{1,0},{1,f[1]}}]}];
line2 = Graphics[{Line[{{6,0},{6,f[6]}}]}];
graph = Plot[f[x],{x,1,6},PlotRange->{{0,6},{0,3.5}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

TrapezoidalRule[f[x],1,6,10]

Compare with Mathematica's NIntegrate procedure.

NIntegrate[f[x],{x,1,6}]

And it can be done with Mathematica's Integrate procedure.

v = Integrate[f[x],{x,1,6}]

v//N

Where the indefinite integral was:

Integrate[f[x],x]

Algorithm 7.2 (Composite Simpson Rule).  To approximate the integral

  aòbÊf(x)dx  ~  h [f(a) + f(b)]/3  +  2h/3  Sk=1k=m-1Êf(x2k)   
                           +  4h/3  Sk=1k=mÊf(x2k-1),

by sampling  f(x)  at the  2m+1  equally spaced points  xk = a + hk  
for  k = 0, 1, 2,...,2m.  Notice that  x0 = a  and  x2m = b.
Section 7.2 Composite Trapezoidal & Simpson's Rule  Page 365

Off[General::spell1];
Clear[SimpsonRule,f0,a,b,m];
SimpsonRule[f0_,a_,b_,m_] :=
Module[{f,k,s,se,so},
  Set@@{f[x_],f0};
  h = (b - a)/(2 m)//N;
  x0 = a//N;
  s = 0;
  se = 0;
  For[k=1,k<=m-1,k++,
    Module[{}, xk = x0 + h 2 k;
               se = se + f[xk]; ]];
  so = 0;
  For[k=1,k<=m,k++,
    Module[{}, xk = x0 + h (2k-1);
               so = so + f[xk]; ]];
  s = h ( f[a] + f[b] + 2 se + 4 so)/3;
  Srule = s;
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to the integral over"];
  Print["[",a,",",b,"] computed by Simpson's rule"];
  Print["using ",2m," subintervals is   ",N[Srule,14]]; ];
On[General::spell1];

Chapter 7, Example 7.6, Page 360.

Integrate  f(x) = 2 + sin( 2 x^1/2 )  over the interval [1,6].

Clear[f,x];
f[x_] = 2 + Sin[2Sqrt[x]]

line1 = Graphics[{Line[{{1,0},{1,f[1]}}]}];
line2 = Graphics[{Line[{{6,0},{6,f[6]}}]}];
graph = Plot[f[x],{x,1,6},PlotRange->{{0,6},{0,3.5}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

SimpsonRule[f[x],1,6,5]

Compare with Mathematica's NIntegrate procedure.

NIntegrate[f[x],{x,1,6}]

And it can be done with Mathematica's Integrate procedure.

v = Integrate[f[x],{x,1,6}]

v//N

Where the indefinite integral was:

Integrate[f[x],x]

Algorithm 7.3  (Recursive Trapezoidal Rule).  To approximate the integral

  aòbÊf(x)dx  ~  h/2  Sk=1k=2^JÊ[f(xk-1) + f(xk)]Ê

by using the trapezoidal rule and successively and increasing number of 
subintervals of  [a,b].  The jth iteration samples f(x) at  2j+1  equally 
spaced points.
Section 7.3 Recursive Rules & Romberg Integration   Page  378

Off[General::spell1];
Clear[SequentialTrapezoid,f0,a,b,n];
SequentialTrapezoid[f0_,a_,b_,n_] :=
Module[{f,i,j,k,s},
  Set@@{f[x_],f0};
  Clear[T];
  T = Table[0,{i,n+1}];
  m = 1;
  h = b - a//N;
  T[[1]] = h (f[a] + f[b])/2//N;
  Print["using 1 \tsubinterval   \t",N[T[[1]],14]];
  For[j=2,j<=n+1,j++,
    Module[{},
      m = 2 m;
      h = h/2;
      s = 0;  
      For[k=1,k<=m/2,k++,
        Module[{}, xk = a + h (2k-1)//N;
                   s = s + f[xk]//N;   ]; ];
      T[[j]] = h s + T[[j-1]]/2; 
      Print["using ",2^(j-1)," \t subintervals  \t",N[T[[j]],14]];]; ];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to the integral over"];
  Print["[",a,",",b,"] computed by sequential trapezoidal rule"];
  Print["using ",2^(j-2)," \t subintervals is   \t",N[T[[n+1]],14]]; ];
On[General::spell1];

Chapter 7, Example 7.11, Page 371.

Integrate  f(x) = 1/x  over the interval [1,5].

Clear[f,x];
f[x_] = 1/x

line1 = Graphics[{Line[{{1,0},{1,f[1]}}]}];
line2 = Graphics[{Line[{{5,0},{5,f[5]}}]}];
graph = Plot[f[x],{x,1,5},PlotRange->{{0,5},{0,1}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

SequentialTrapezoid[f[x],1,5,10]

Compare with Mathematica's integrals.

v = Integrate[f[x],{x,1,5}]

N[v,14]

Where the indefinite integral was:

Integrate[f[x],x]

Algorithm 7.4   (Romberg Integration).  To approximate the integral

  aòbÊf(x)dx  ~  R(j,j) 

by generating a table of approximations R(i,k) for  k <= i, and using R(j,j) as the final answer.  The approximations R(j,k) are stored in a special lower triangular matrix.  The elements R(j, 0) of column 0 are computed using the sequential Trapezoidal rule based on  2j  subintervals of  [a,b];  then R(j,k) is computed using  Romberg's rule.  Elements of row  j  are

R(j,k)  =  R(j,k-1)  +  (R(j,k-1)-R(j-1,k-1))/(4kÊ-Ê1)Ê    for  1 <= k <=  j. 

The algorithm is terminated in the j-th row when
|R(j,j) - R(j-1,j-1)| < Tol.
Section	7.3	Recursive Rules & Romberg Integration   Page  379

At least four rows of the table are computed. The algorithm is terminated in the jth row when the approximations  R(j-1,j-1) and R(j,j) differ by less than the preassigned value  Tol.

Off[General::spell1];
Clear[Romberg,f0,a,b,Tol,digits];
Romberg[f0_,a_,b_,Tol_,digits_] :=
Module[{f,k},
  Set@@{f[x_],f0};
  SeqTrapRule[j0_] :=
  Module[{J,p},    
    h = h/2;
    J = j0;
    sum = 0;
    For[p=1,p<=m,p++,
      Module[{},
        xp = a + h (2p-1)//N;
        sum = sum + f[xp]//N; ]; ];
    R[[J+1,0+1]] = h sum + R[[J-1+1,0+1]]/2//N;
    R[[J+1,0+1]] = N[R[[J+1,0+1]],digits];
    m = 2m; ];
  min = 1;
  max = 8;
  m = 1;
  close = 1;
  j = 0;
  h = b - a//N;
  Clear[R,ii,jj];
  R = Table[0,{ii,15},{jj,ii}];
  R[[0+1,0+1]] = h (f[a] + f[b])/2//N;
  Print[N[R[[0+1,0+1]],14]];
  R[[0+1,0+1]] = N[R[[0+1,0+1]],digits];
  While[((j<max)) || (j<min),
    j = j+1;
    SeqTrapRule[j];
    For[k=1,k<=j,k++,
      Block[{},
      R[[j+1,k+1]] = R[[j+1,k-1+1]] + 
        (R[[j+1,k-1+1]] - R[[j-1+1,k-1+1]])/(4^k-1);
      R[[j+1,k+1]] = N[R[[j+1,k+1]],digits];   ]; ];
    Print[N[R[[j+1,j+1]],14]];
    close = Abs[R[[j-1+1,j-1+1]] - R[[j+1,j+1]]];  ];
  Rrule = R[[j+1,j+1]];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to the integral over"];
  Print["[",a,",",b,"] computed with Romberg integration"];
  Print["using ",2^j," subintervals is   ",N[Rrule,14]]; ];
On[General::spell1];

Chapter 7, Example 7.14, Page 377.

Integrate  f(x) = (x^2 + x + 1)cos(x)  over the interval [0,p/2].

Clear[f,x];
f[x_] = (x^2 + x + 1) Cos[x]

line1 = Graphics[{Line[{{0,0},{0,f[0]}}]}];
line2 = Graphics[{Line[{{Pi/2,0},{Pi/2,f[Pi/2]}}]}];
graph = Plot[f[x],{x,0,Pi/2},PlotRange->{{0,1.75},{0,1.75}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

Romberg[f[x],0,Pi/2,10^-8,13]

Compare with Mathematica's integrals.

v = Integrate[f[x],{x,0,Pi/2}]

v = N[v,14]

Where the indefinite integral was:

Integrate[f[x],x]

Algorithm 7.5  (Adaptive Quadrature Using Simpson's Rule). 
To approximate the integral

aòbÊf(x)dx  ~
   Sk=1k=mÊhk/3 [f(x4k-4) + 4 f(x4k-3) + 2f(x4k-2) + 4 f(x4k-1) + f(x4k)].

The composite Simpson rule is applied to the 4m subintervals 
[x4k-4,x4k]  where [a,b] = [x0,x4m]  and  x4k-4+j = x4k-4 + j hk  
for each  k = 1,  2, ...,m  and  j = 1,..4. 
Section	7.4	Adaptive Quadrature   Page  389

Adaptive Quadrature: Simpson's Rule

First, execute the following group of cells:

Off[General::spell1];
Clear[SimpRule,a0,b0,tol0];
SimpRule[a0_,b0_,tol0_] :=
  Module[{a,c,b,fa,fc,fb,s,t,err,tol,tocheck},
  a = a0 //N;
  b = b0 //N;
  h = (b - a)/2 //N;
  c = a + h;
  fa = f[a] //N;
  fc = f[c] //N;
  fb = f[b] //N;
  s = h (fa + 4fc + fb)/3;
  t = s;
  tol = tol0;  
  err = tol0;  
  Return[{a,c,b,fa,fc,fb,s,t,err,tol,tocheck}]; ];

Clear[Refine,SR0,p];
Refine[SR0_,p_] :=
  Block[{a,c,b,fa,fc,fb,s,t,err,tol,tocheck},
  state = done;
  SR  = SR0;
  sr0 = SR0[[p]];
  {a,c,b,fa,fc,fb,s,t,err,tol,tocheck} = sr0;
  If[tocheck==pass, Return[SR]];
  sr1 = SimpRule[a,c,tol/2];
  sr2 = SimpRule[c,b,tol/2];
  err = Abs[sr0[[7]]-sr1[[7]]-sr2[[7]]]/10;
  If[  err<tol, sr0[[11]] = pass];
  If[  err<tol, SR[[p]] = sr0;
    SR[[p,8]] = sr1[[7]] + sr2[[7]];
    SR[[p,9]] = err,
    SR[[p]] = sr2;
    SR = Insert[SR,sr1,p];
    state = iterating ];
  Return[SR]; ];

Clear[RefineBoole,SR0,p];
RefineBoole[SR0_,p_] :=
  Block[{a,c,b,fa,fc,fb,s,t,err,tol,tocheck},
  state = done;
  SR  = SR0;
  sr0 = SR0[[p]];
  {a,c,b,fa,fc,fb,s,t,err,tol,tocheck} = sr0;
  If[tocheck==pass, Return[SR]];
  sr1 = SimpRule[a,c,tol/2];
  sr2 = SimpRule[c,b,tol/2];
  err = Abs[sr0[[7]]-sr1[[7]]-sr2[[7]]]/10;
  If[  err<tol, sr0[[11]] = pass];
  If[  err<tol, SR[[p]] = sr0;
    SR[[p,8]] = (16(sr1[[7]] + sr2[[7]])-sr0[[7]])/15;
    SR[[p,9]] = err,
    SR[[p]] = sr2;
    SR = Insert[SR,sr1,p];
    state = iterating ];
  Return[SR]; ];

Clear[AdaptQuad,f0,a0,b0,tol0];
AdaptQuad[f0_,a0_,b0_,tol0_] := Block[{},
  f[x_] = f0;
  SR = {SimpRule[a0,b0,tol0]};
  state = iterating;
  While[state==iterating,
    For[i=Length[SR],1<=i,i--, Refine[SR,i]] ];
  RS = Transpose[SR];
  RSxsy = {RS[[1]],RS[[3]],RS[[8]],
           RS[[9]],RS[[10]]};
  SRxsy = Transpose[RSxsy];
  YX = {RS[[3]],RS[[8]]};
  XY = Transpose[YX];
  XY = Prepend[XY,{RS[[1,1]],0}];
  For[i=2,i<=Length[XY],i++,
    XY[[i,2]] = XY[[i-1,2]] + XY[[i,2]] ];
  Clear[i];
  sum = Sum[SRxsy[[i,3]],{i,1,Length[SR]}];
  Return[N[sum,12]] ];

Clear[AdaptBoole,f0,a0,b0,tol0];
AdaptBoole[f0_,a0_,b0_,tol0_] := Block[{},
  f[x_] = f0;
  SR = {SimpRule[a0,b0,tol0]};
  state = iterating;
  While[state==iterating,
    For[i=Length[SR],1<=i,i--, RefineBoole[SR,i]] ];
  RS = Transpose[SR];
  RSxsy = {RS[[1]],RS[[3]],RS[[8]],
           RS[[9]],RS[[10]]};
  SRxsy = Transpose[RSxsy];
  YX = {RS[[3]],RS[[8]]};
  XY = Transpose[YX];
  XY = Prepend[XY,{RS[[1,1]],0}];
  For[i=2,i<=Length[XY],i++,
    XY[[i,2]] = XY[[i-1,2]] + XY[[i,2]] ];
  Clear[i];
  sum = Sum[SRxsy[[i,3]],{i,1,Length[SR]}];
  Return[N[sum,12]] ];
On[General::spell1];

Chapter 7, Exercise 7.16, Page 388.

Integrate  f(x) = 13(x - x^2)exp(-3x/2)  over the interval [0,4].

Clear[f,x];
f[x_] = 13(x-x^2)Exp[-3x/2]

line1 = Graphics[{Line[{{0,0},{0,f[0]}}]}];
line2 = Graphics[{Line[{{4,0},{4,f[4]}}]}];
graph = Plot[f[x],{x,0,4},PlotRange->{{0,4},{-1.4,1.8}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

ans1 = AdaptQuad[f[x],0,4,0.00001]

ans2 = AdaptBoole[f[x],0,4,0.00001]

Use Mathematica's NIntegrate procedure.

N[NIntegrate[f[x],{x,0,4}],20]

Compare with Mathematica's integrals.

v = Integrate[f[x],{x,0,Pi/2}]

v = N[v,14]

Where the indefinite integral was:

Integrate[f[x],x]

Algorithm 7.6  (Gauss-Legendre Quadrature).  
To approximate the integral

aòbÊf(x)dx  ~  (b-a)/2  [wn,1f(tn,1) + wn,2f(tn,2) + ... + wn,nf(tn,n)]

by sampling f(x) at the n unequally spaced points tn,1, tn,2,..., tn,n. The 
change of variable t = (a+b)/2 + x(b-a)/2 and dt = (b-a)/2 dx  is used.  
The abscissas xn,1, xn,2, ..., xn,n  and the corresponding weights 
wn,1, wn,2, ..., wn,n  must be obtained from a table of known values.
Section	7.5	Gauss-Legendre Integration   Page  397

First, execute the following group of cells:

Off[General::spell1];
Clear[GaussQuad,f0,a0,b0,Tol,digits];
GaussQuad[f0_,a0_,b0_,Tol_,digits_] :=
Module[{close,f,j,k,mid,xx,sum,wide},
  Set@@{f[x_],f0};
  mid  = (a0 + b0)/2;
  wide = (b0 - a0)/2;
  close = 1;
  j = 1;
  xx = mid + a[1,1];
  Clear[i,Q];
  Q = Table[0,{i,8}];  
  Q[[1]] = w[1,1] f[xx] wide;
  Print[N[Q[[1]],14]];
  While[((close>Tol) || (j<3)) && (j<8),
    j = j+1;
    sum = 0;
    For[k=1,k<=j,k++,
      Module[{xx},
        xx = mid + a[j,k] wide;
        sum =sum + w[j,k] f[xx]; ]; ];
    Q[[j]] = sum wide;
    Print[N[Q[[j]],14]];
    close  = Abs[Q[[j]] - Q[[j-1]]];  ]; 
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to the integral over"];
  Print["[",a0,",",b0,"] computed with Gauss Legendre integration"];
  Print["using ",j," sample points is   ",N[Q[[j]],14]]; ];

SetUpWeights :=
Block[{j,k},
  For[j=1,j<=8,j++,
    For[k=1,k<=Floor[j/2],k++,
      Block[{},
        a[j,j+1-k] = -a[j,k];
        w[j,j+1-k] = w[j,k]; ]; ]; ]; ];

Off[General::spell1];

Clear[a,w];
a[1,1] = 0.0;
w[1,1] = 2.0;
a[2,1] = 0.577350269189626;
w[2,1] = 1.0;
a[3,1] = 0.774596669241483;
a[3,2] = 0.0;
w[3,1] = 0.555555555555556;
w[3,2] = 0.888888888888889;
a[4,1] = 0.861136311594053;
a[4,2] = 0.339981043584856;
w[4,1] = 0.347854845137454;
w[4,2] = 0.652145154862546;
a[5,1] = 0.906179845938664;
a[5,2] = 0.538469310105683;
a[5,3] = 0.0;
w[5,1] = 0.236926885056189;
w[5,2] = 0.478628670499366;
w[5,3] = 0.568888888888889;
a[6,1] = 0.932469514203152;
a[6,2] = 0.661209386466265;
a[6,3] = 0.238619186083197;
w[6,1] = 0.171324492379170;
w[6,2] = 0.360761573048139;
w[6,3] = 0.467913934572691;
a[7,1] = 0.949107912342759;
a[7,2] = 0.741531185599394;
a[7,3] = 0.405845151377397;
a[7,4] = 0.0;
w[7,1] = 0.129484966168870;
w[7,2] = 0.279705391489277;
w[7,3] = 0.381830050505119;
w[7,4] = 0.417959183673469;
a[8,1] = 0.960289856497536;
a[8,2] = 0.796666477413627;
a[8,3] = 0.525532409916329;
a[8,4] = 0.183434642495650;
w[8,1] = 0.101228536290376;
w[8,2] = 0.222381034453374;
w[8,3] = 0.313706645877887;
w[8,4] = 0.362683783378362;

Chapter 7, Exercise 3, Page 398.

Integrate  f(x) = sin(x)/x  over the interval [0,1].

First you must set up the weights and enter the function.

SetUpWeights
Clear[f,x];
f[x_] = Sin[x]/x

line1 = Graphics[{Line[{{0.00001,0},{0.00001,f[0.00001]}}]}];
line2 = Graphics[{Line[{{1,0},{1,f[1]}}]}];
graph = Plot[f[x],{x,0.00001,1},PlotRange->{{0,1},{0,1}},
                          DisplayFunction->Identity];
Show[graph,line1,line2,
     DisplayFunction->$DisplayFunction];

GaussQuad[f[x],0,1,0.00001,19]

Compare with Mathematica's integrals.

v = Integrate[f[x],{x,0,1}]

v = N[v,14]

Where the indefinite integral was:

Integrate[f[x],x]