Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

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.

<<  Graphics`FilledPlot`Clear["Global`*"];
Off[General::"spell1"];
Clear[a, a0, AdaptBoole, AdaptQuad, b, b0, c, close, digits, err, f, 
        fa, fb, fc, f0, graph, gr1, GaussQuad, i, j, J, k, line1, line1, m, 
        mid, n, p, Refine, RefineBoole, Romberg, s, se, so, sum, 
        SequentialTrapezoid, SimpRule, SimpsonRule, SR0, t, tol, tol0, 
        tocheck, Tol, TrapezoidalRule, w, wide, x, xx];
On[General::"spell1"];

Algorithm 7.1  Composite Trapezoidal Rule. To approximate the 
	(∫_a^b)f (x) dx  ≈  (h/2)[f (a) + f (b)] + h  (∑_(k = 1)^(m - 1))f ((Subscript[x, k]))
by sampling  f (x)  at the  m+1 equally spaced points  (Subscript[x, k]) = a + k h  
for  k = 0, 1, ⋯, m.   Notice that  (Subscript[x, 0]) = a  and  (Subscript[x, m]) = b.
Section  7.2   Composite Trapezoidal & Simpson's Rule   Page 365

Execute this cell to activate Composite Trapezoidal Rule.
TrapezoidalRule[f0_,a_,b_,m_]
f0_ is the function f(x),
a_ is the left endpoint,
b_ is the right endpoint,
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[TrapezoidalRule, f0, a, b, m]; \(TrapezoidalRule[f0_, a_, b_, m_] := 
	
      Module[{f, k}, 
		Set@@{f[x_], f0}; 
		h = N[\(b - a\)\/m]; 
		
        x0 = N[a]; 
		s = 0; 
		
        For[k = 1, k ≤ m - 1, \(k++\), 
			
          Module[{}, 
				xk = x0 + h\ k; 
				
            s = N[s + f[xk]]; ]]; 
		
        s = N[1\/2\ h\ \((f[a] + f[b])\) + h\ s]; 
		Trule = s; 
		
        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 (Sqrt[x]))  over the interval  [1, 6].

\(Clear[f, graph, gr1, line1, line2, x]; 
f[x_] = 2 + Sin[2\ \@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.0}}, 
        DisplayFunction -> Identity]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 1, 6}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[graph, gr1, line1, line2, DisplayFunction -> $DisplayFunction]; 
    
Print["The area under y = ", f[x]]; \)

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

Compare with Mathematica's NIntegrate procedure.

\(nval\  = \ NIntegrate[f[x], {x, 1, 6}]; 

    Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x; = "", 
      N[nval, 14]]; \)

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

\(val\  = \ ∫\_1\%6 f[x] \[DifferentialD]x; 

    Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x = "", val]; 
    
Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x = "", 
      N[val, 14]]; \)

Where Mathematica used the indefinite integral:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];

Algorithm 7.2 Composite Simpson Rule.  To approximate the integral  (∫_a^b)f (x) dx ≈ (h/3)[f (a) + f (b)]+((2 h)/3) (∑_(k = 1)^(m - 1))f ((Subscript[x, 2  k]))+((4 h)/3)(∑_(k = 1)^m)f ((Subscript[x, 2  k - 1]))
by sampling  f (x)  at the  2  m + 1  equally spaced points  (Subscript[x, k]) = a + k h for  k = 0, 1, ⋯, 2  m.  Notice that  (Subscript[x, 0]) = a  and  (Subscript[x, 2  m]) = b.
Section  7.2   Composite Trapezoidal & Simpson's Rule  Page 365

Execute this cell to activate Composite Simpson Rule.
SimpsonRule[f0_,a_,b_,m_]
f0_ is the function f(x),
a_ is the left endpoint,
b_ is the right endpoint,
m_  here 2m is the number of subintervals.

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 = N[\(b - a\)\/\(2\ m\)]; 
		x0 = N[a]; 
		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\ \((2\ k - 1)\); 
				
            so = so + f[xk]; ]]; 
		
        s = 1\/3\ h\ \((f[a] + f[b] + 2\ se + 4\ so)\); 
		Srule = s; 

        		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 ", 2\ m, " subintervals is   ", 
          N[Srule, 14]]; ]; 
On[General::"spell1"]; \)

Chapter 7, Example 7.6, Page 360.
Integrate  f (x)  =  2 + sin (2 (Sqrt[x]))  over the interval  [1, 6].

\(Clear[f, graph, gr1, line1, line2, x]; 
f[x_] = 2 + Sin[2\ \@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.0}}, 
        DisplayFunction -> Identity]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 1, 6}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[graph, gr1, line1, line2, DisplayFunction -> $DisplayFunction]; 
    
Print["The area under y = ", f[x]]; \)

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

Compare with Mathematica's NIntegrate procedure.

\(nval\  = \ NIntegrate[f[x], {x, 1, 6}]; 

    Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x; = "", 
      N[nval, 14]]; \)

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

\(val\  = \ ∫\_1\%6 f[x] \[DifferentialD]x; 

    Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x = "", val]; 
    
Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x = "", 
      N[val, 14]]; \)

Where Mathematica used the indefinite integral:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];

Algorithm 7.3  Recursive Trapezoidal Rule.  To approximate the integral
	(∫_a^b)f (x) dx  ≈  (h/2) (∑_(k = 1)^((2^j)))[f ((Subscript[x, k - 1]))  + f ((Subscript[x, k]))]
by using the trapezoidal rule and successively and increasing number of subintervals 
of  [a, b].  The  (j^th)  iteration samples  f (x)  at  (2^j) + 1  equally spaced points.
Section  7.3   Recursive Rules & Romberg Integration   Page  378

Execute this cell to activate Recursive Trapezoidal Rule.
SequentialTrapezoid[f0_,a_,b_,n_]
f0_ is the function f(x),
a_ is the left endpoint,
b_ is the right endpoint,
n_  is the number of recursions.

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 = N[b - a]; 
		
        T[[1]] = 
          N[1\/2\ h\ \((f[a] + f[b])\)]; 
		
        Print["using 1 	subinterval get     	", 
          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 = N[a + h\ \((2\ k - 1)\)]; 

                						s = N[s + f[xk]]; ]; \)]; 
				
            T[[j]] = 
              h\ s + 1\/2\ T[[j - 1]]; 

            				Print["Using ", 2\^\(j - 1\), 
              " 	 subintervals get  	", 
              N[T[[j]], 14]]; ]; \)]; 

        		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\), " 	 subintervals is   	", 
          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, graph, gr1, line1, line2, 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]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 1, 5}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[graph, gr1, line1, line2, DisplayFunction -> $DisplayFunction]; 
    
Print["The area under y = ", f[x]]; \)

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

Compare with Mathematica's NIntegrate procedure.

\(nval\  = \ NIntegrate[f[x], {x, 1, 5}]; 

    Print[\*"" \!\(∫\_1\%6\)f[x]\[DifferentialD]x; = "", 
      N[nval, 14]]; \)

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

\(val\  = \ ∫\_1\%5 f[x] \[DifferentialD]x; 

    Print[\*"" \!\(∫\_1\%5\)f[x]\[DifferentialD]x = "", val]; 
    
Print[\*"" \!\(∫\_1\%5\)f[x]\[DifferentialD]x = "", 
      N[val, 14]]; \)

Where Mathematica used the indefinite integral:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];

Algorithm 7.4   Romberg Integration.  To approximate theintegral
		(∫_a^b)f (x) dx  ≈ (Subscript[R, j, j])
by generating a table of approximations  (Subscript[R, i, k])  for  k ≤ i, and using  (Subscript[R, j, j])  as the final answer.  The approximations  (Subscript[R, j, k])  are stored in a special lower triangular matrix.  
The elements  (Subscript[R, j, 0])  of column 0 are computed using the sequential Trapezoidal rule 
based on  (2^j)  subintervals of  [a, b];  then  (Subscript[R, j, k])  is computed using Romberg's rule.  Elements of row  j  are
	(Subscript[R, j, k]) = (Subscript[R, j, k - 1]) + (((Subscript[R, j, k - 1])   -   (Subscript[R, j - 1, k - 1]))/((4^k)  -  1))   for  1≤ k ≤  j. 
The algorithm is terminated in the j-th row when  (Subscript[(| R), j, j])  -  (Subscript[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 row  j  when the approximations   (Subscript[R, j - 1, j - 1])  and  (Subscript[R, j, j])  differ by less than the preassigned value  Tol.

Execute this cell to activate Romberg Integration.
Romberg[f0_,a_,b_,Tol_,digits_]
f0_ is the function f(x),
a_ is the left endpoint,
b_ is the right endpoint,
Tol_ is the tolerance,
digits_ is the number of digits displayed.

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 = N[a + h\ \((2\ p - 1)\)]; 

                						sum = N[sum + f[xp]]; ]; \)]; 
				
            R[[J + 1, 0 + 1]] = 
              N[h\ sum + 
                  1\/2\ R[[J - 1 + 1, 
                      0 + 1]]]; 
				
            R[[J + 1, 0 + 1]] = 
              N[R[[J + 1, 0 + 1]], 
                digits]; 
				m = 2\ m; ]; 
		min = 1; 
		max = 8; 
        
		m = 1; 
		close = 1; 
		j = 0; 
		h = N[b - a]; 
		
        Clear[R, ii, jj]; 
		R = Table[0, {ii, 15}, {jj, ii}]; 
		
        R[[0 + 1, 0 + 1]] = 
          N[1\/2\ h\ \((f[a] + f[b])\)]; 
		
        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,(π/2)].

\(Clear[f, graph, gr1, line1, line2, x]; 

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

    line1 = Graphics[{Line[{{0, 0}, {0, f[0]}}]}]; 

    \(line2 = Graphics[{Line[{{π\/2, 0}, {π\/2, f[π\/2]}}]}]; 
    \)\)\(graph = 
      Plot[f[x], {x, 0, π\/2}, PlotRange -> {{0, 1.75}, {0, 1.75}}, 
        DisplayFunction -> Identity]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 0, π\/2}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[graph, gr1, line1, line2, DisplayFunction -> $DisplayFunction]; 
    
Print["The area under y = ", f[x]]; \)

Romberg[f[x], 0, (π/2), 0.00000001, 13];

Compare with Mathematica's NIntegrate procedure.

\(nval\  = \ NIntegrate[f[x], {x, 0, π\/2}]; 

    Print[\*"" (∫_0^((π/2)))f[x]\[DifferentialD]x; = "", N[nval, 14]]; \)

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

\(val\  = \ ∫\_0\%\(π\/2\)f[x] \[DifferentialD]x; 

    Print[\*"" (∫_0^((π/2)))f[x]\[DifferentialD]x = "", val]; 

    Print[\*"" (∫_0^((π/2)))f[x]\[DifferentialD]x = "", N[val, 14]]; \)

Where Mathematica used the indefinite integral:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];

Algorithm 7.5  Adaptive Quadrature Using Simpson's Rule. 
To approximate the integral
(∫_a^b)f (x) dx ≈
(h/3) (∑_(k = 1)^m)[f ((Subscript[x, 4  k - 4])) + 4  f ((Subscript[x, 4  k - 3])) + 
                2  f ((Subscript[x, 4  k - 2])) + 4  f ((Subscript[x, 4  k - 1])) + 
                f ((Subscript[x, 4  k]))]
The composite Simpson rule is applied to the  4m  subintervals 
[(Subscript[x, 4  k - 4]), (Subscript[x, 4  k])]  where  [a, b] = [(Subscript[x, 0]), (Subscript[x, 4  m])]  and  (Subscript[x, 4  k - 4 + j]) = (Subscript[x, 4  k - 4]) + j (Subscript[h, k])
for each  k = 1, 2, ⋯, m  and  j = 1, ⋯, 4. 
Section  7.4   Adaptive Quadrature   Page  389

Execute this cell to activate Adaptive Quadrature Using Simpson's Rule.
AdaptQuad[f0_,a0_,b0_,tol0_]
AdaptBoole[f0_,a0_,b0_,tol0_]
f0_ is the function f(x),
a0_ is the left endpoint,
b0_ is the right endpoint,
tol0_ is the tolerance.

\(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\  + \ 4  fc\  + \ 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[{}, 
\ \ 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[{}, 
\ \ 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)) (e^(((-3 x)/2)))  over the interval  [0, 4].

\(Clear[f, graph, gr1, gr2, line1, line2, x]; 

    f[x_] = 13\ \((x - x\^2)\)\ E\^\(\(\(-3\)\ x\)\/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]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 0, 1}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    gr2\  = \ FilledPlot[{f[x]}, 
\ {x, \ 1, 4}, 
        Fills -> {{{1, Axis}, RGBColor[1, 0, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[gr1, gr2, graph, line1, line2, 
      DisplayFunction -> $DisplayFunction]; 

    Print["The area under y = ", f[x]]; \)

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 = (∫_0^4) f[x] \[DifferentialD]x

v = N[v, 14]

Where the indefinite integral was:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];

Algorithm 7.6  Gauss-Legendre Quadrature.  To approximate the integral
	(∫_a^b)f (t) dt  ≈  ((b - a)/2) (∑_(k = 1)^n)(Subscript[w, n, k]) f ((Subscript[t, n, k]))
by sampling  f (x)  at the  n  unequally spaced points  (Subscript[t, n, 1]), (Subscript[t, n, 2]), ⋯, (Subscript[t, n, n]).
The change of variable t=((a + b)/2)+((b - a)/2)x  and  dt=((b - a)/2)dx  is used.  
The abscissas  (Subscript[x, n, 1]), (Subscript[x, n, 2]), ⋯, (Subscript[x, n, n])  and the corresponding weights 
(Subscript[w, n, 1]), (Subscript[w, n, 2]), ⋯, (Subscript[w, n, n])  must be obtained from a table of known values.
Section  7.5   Gauss-Legendre Integration   Page  397

Execute this cell to activate Gauss-Legendre Quadrature.
GaussQuad[f0_,a0_,b0_,Tol_,digits_]
f0_ is the function f(x),
a0_ is the left endpoint,
b0_ is the right endpoint,
Tol_ is the tolerance,
digits_ is the number of digits displayed.

Off[General::"spell1"]; \(Clear[a, GaussQuad, f0, a0, b0, Tol, digits, w]; 

    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]; ]; \)]; \)]; 

        Print["The weights and absicssas have been entered."]; ]; \)\)

Off[General::"spell1"];

Execute this cell to activate The weights and abscissas.

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, graph, gr1, line1, line2, 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.1}, {0, 1.1}}, 
        DisplayFunction -> Identity]; 

    gr1\  = \ FilledPlot[{f[x]}, 
\ {x, 0.00001, 1}, 
        Fills -> {{{1, Axis}, RGBColor[0, 1, 0]}}, 
        DisplayFunction -> Identity]; 

    Show[graph, gr1, line1, line2, DisplayFunction -> $DisplayFunction]; 
    
Print["The area under y = ", f[x]]; \)

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

Compare with Mathematica's integrals.

v = (∫_0^1) f[x] \[DifferentialD]x

v = N[v, 14]

Where the indefinite integral was:

Print[" ∫f[x]\[DifferentialD]x= ", 
      ∫f[x] \[DifferentialD]x];