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 6.   Numerical Differentiation
	Algorithm 6.1  Differentiation Using Limits.
	Algorithm 6.2  Differentiation Using Extrapolation.
	Algorithm 6.3  Differentiation Based on n+1 Nodes.

Clear["Global`*"];
Off[General::"spell1"];
Clear[big, d, dot, dots, DerAtX1, e, Epsilon, ExtrapolateDQ, f, f0, 
        h0, j, k, line, LimitDQ, max, min, n, small, Tol, x0, x, xy, X, XY];
On[General::"spell1"];

Algorithm 6.1 Differentiation Using Limits.  
To approximate  f'(x)  numerically by generating the sequence 
     f'(x) ≈ (Subscript[(d ), k])= ((f (x + (h/(2^k))) - f (x - (h/(2^k))))/((h/(2^(k - 1))))),   for  k = 0,1,...,n   
until  (Subscript[(| d ), n + 1]) -  (Subscript[(d ), n]) 
        |  
         ≤  
          |(Subscript[( d ), n]) -  (Subscript[(d ), n - 1]) 
            |   or  (Subscript[(| d ), n]) -  (Subscript[(d ), n - 1]) 
        |   ≺   Tol,   
which is an attempt to find the best  approximation   f'(x) ≈ (Subscript[(d ), n]) .
Section  6.1   Approximating the Derivative   Page  326

Execute this cell to activate Differentiation Using Limits.
LimitDQ[f0_,x0_,h0_,Tol_]
f0_ is the function f(x),
x0_ is the fixed value of x,
h0_ is the step size,
Tol_ is the tolerance.

Off[General::"spell1"]; Clear[LimitDQ, f0, x0, h0, Tol]; LimitDQ[f0_, x0_, h0_, Tol_]:=
	Module[{f, j, big, max, min, small},
		Set@@{f[x_], f0};
		small = N[(1/(10^7))];
		big = N[(10^6)];
		max = 32;
		min = 2;
		Clear[d, e];
		d = {};
		e = {1};
		h = N[h0];
		dydx = N[((0.5 (f[x0 + h] - f[x0 - h]))/h)];
		d = Append[d, dydx];
		n = 2;
		\(For[j = 2, j ≤ 3, \(j++\), 
			h = h\/2; 

                			dydx = N[\(0.5\ \((f[x0 + h] - f[x0 - h])\)\)\/h]; 

                			Print[N[dydx, 14], "	     h = ", h]; 
			
                d = Append[d, dydx]; 
			n = Length[d]; 
			
                Err = Abs[
                    d[[n]] - 
                      d[[n - 1]]]; 

                			e = Append[e, Err]; 
			
                RelErr = 
                  \(2\ Err
                      \)\/\(Abs[d[[n]]]
                         + Abs[
                        d[[n - 1]]] + 
                      small\)]\);
		\(While[\((
                      e[[n]] > 
                          e[[n - 1]] || 
                        RelErr > Tol)\) && n < max && 
                    Abs[d[[n]]] < big || 
                  n < min, 
			h = h\/2; 
			
                dydx = N[\(0.5\ \((f[x0 + h] - f[x0 - h])\)\)\/h]; 
			
                Print[N[dydx, 14], "	     h = ", h]; 
			
                d = Append[d, dydx]; 
			n = Length[d]; 
			
                Err = Abs[
                    d[[n]] - 
                      d[[n - 1]]]; 

                			e = Append[e, Err]; 
			
                RelErr = 
                  \(2\ Err
                      \)\/\(Abs[d[[n]]]
                         + 
                      Abs[d[[n - 1]]] + 
                      small\)]\);
		Print[" "];
		Print["The function is f[x] = ", f[x]];
		Print["A numerical approximation to 
f'[",x0,"] is:"];
		Print["
f'[",x0,"] ≈ ",N[d[[n]], 14]
                    ];];
On[General::"spell1"];

Chapter 6, Example 6.2, Page 320.
Use formula (3) on page 318 to approximate 
the derivative of  f(x) = cos(x)  at  x = 1.

\(f[x_] = Cos[x]; 
LimitDQ[f[x], 1, 1, 0.0000001]; \)

x0  = 1;
y0  = f[x0];
d0  =  d[[n]];
grf = Plot[f[x], {x, 0.0, 1.5}, 
            PlotStyle -> RGBColor[0, 1, 0], 
            DisplayFunction -> Identity];
dot = ListPlot[{{x0, y0}}, PlotStyle -> {PointSize[0.025]}, 
            DisplayFunction -> Identity];
\(line = 
          Graphics[{RGBColor[1, 0, 0], 
			
              Line[{{x0 - 1\/2, y0 - d0\/2}, {x0 + 1\/2, y0 + d0\/2}}]}]\);
Show[grf, dot, line, PlotRange -> {{0, 1.5}, {0, 1}}, 
          DisplayFunction -> $DisplayFunction, 
          AxesLabel -> {"x", "y    "}];
Print["The function is f[x] = ", f[x]];
Print["
f'[",x0,"] ≈ ",N[d[[n]], 14]];

Compare with the analytic method of computing derivatives.

Print["
f'[1] = ",N[d[[n]], 14]];
Print["
f'[1] = ",N[f'[1],14]];

Algorithm  6.2  Differentiation Using Extrapolation.  
To approximate  f'(x) numerically by generating a table of approximations  d (j, k)  for  k ≤ j,  and using  f'(x) ≈ d (n, n)  as the final answer.  The approximations  d (j, k)  are stored in a lower triangular matrix.  The first row is       
          d (j, 0) = ((f (x + (h/(2^j))) - f (x - (h/(2^j))))/((h/(2^(j - 1)))))
and the elements in row  j  are:
d (j, k) = d (j, k - 1) + ((d (j, k - 1) - d (j - 1, k - 1))/((4^k)  -  1)),  for 1 ≤ k ≤  j.
Section  6.1   Approximating the Derivative   Page  327

Execute this cell to activate Differentiation Using Extrapolation.
ExtrapolateDQ[f0_,x0_,h0_,Tol_,Epsilon_]
f0_ is the function f(x),
x0_ is the fixed value for x,
h0_ is the step size,
Tol_ is the relative error tolerance,
Epsilon_ is the error bound.

Off[General::"spell1"]; Clear[ExtrapolateDQ, f0, x0, h0, Tol, Epsilon]; ExtrapolateDQ[f0_, x0_, h0_, Tol_, Epsilon_]:=
	Module[{f, j, big, max, min, small},
		Set@@{f[x_], f0};
		small = N[(1/(10^7))];
		big = N[(10^6)];
		max = 32;
		min = 2;
		Clear[d, e];
		d = Table[0, {11}, {11}];
		h = N[h0];
		j = 1;
		Error = 1.0;
		RelErr = 1.0;
		d[[1, 1]] = 
                N[((0.5 (f[x0 + h] - f[x0 - h]))/h)];
		Print[
                N[d[[1, 1]], 14], 
                "\t     h = ", h];
		\(While[RelErr > Tol && Error > Epsilon && j < 10, h = h\/2; 

                			d[[j + 1, 1]] = 
                  N[\(0.5\ \((f[x0 + h] - f[x0 - h])\)\)\/h]; 
			
                For[k = 1, k ≤ j, \(k++\), 
				
                  \(d[[j + 1, k + 1]] = 
                    d[[j + 1, 
                        k - 1 + 1]] + 
                      \(d[[j + 1, k]] - 
                          d[[j, k
                            ]]\)\/\(4\^k - 1\); \)]; 

                			Print[
                  N[d[[j + 1, j + 1]], 
                    14], "	     h = ", h]; 
			
                Error = 
                  N[Abs[d[[j + 1, 
                          j + 1]] - 
                        d[[j, j]]]]; 

                			RelErr = 
                  N[\(2\ Error
                        \)\/\(Abs[
                          d[[j + 1, 
                            j + 1]]] + 
                        Abs[d[[j, j]]] + 
                        small\)]; 
			j = j + 1; 
			n = j; ]\);
		Print[" "];
		Print["The function is f[x] = ", f[x]];
		Print["A numerical approximation to 
f'[",x0,"] is:"];
		Print["
f'[",x0,"] ≈ ",N[d[[n, n]], 14]
                    ];];
On[General::"spell1"];

Chapter 6, Example 6.3, Page 325.
Use extrapolation to approximate the derivative of  f(x) = cos(x)  at  x = 1.

\(f[x_] = Cos[x]; 
ExtrapolateDQ[f[x], 1, 1, 0.0000001, 0.0000001]; \)

x0  = 1;
y0  = f[x0];
d0  =  d[[n, n]];
grf = 
          Plot[f[x], {x, 0.0, 1.5}, PlotStyle -> RGBColor[0, 1, 0], 
            DisplayFunction -> Identity];
dot = ListPlot[{{x0, y0}}, PlotStyle -> {PointSize[0.025]}, 
            DisplayFunction -> Identity];
\(line = 
          Graphics[{RGBColor[1, 0, 0], 
			
              Line[{{x0 - 1\/2, y0 - d0\/2}, {x0 + 1\/2, y0 + d0\/2}}]}]\);
Show[grf, dot, line, PlotRange -> {{0, 1.5}, {0, 1}}, 
          DisplayFunction -> $DisplayFunction, 
          AxesLabel -> {"x", "y    "}];
Print["The function is f[x] = ", f[x]];
Print["
f'[",x0,"] ≈ ",N[d[[n, n]], 14]];

Compare with the analytic method of computing derivatives.

Print["
f'[1] = ",N[d[[n, n]], 14]];
Print["
f'[1] = ",N[f'[1],14]];

Algorithm 6.3  Differentiation Based on n+1 Nodes. 
To approximate  f'(x)  numerically by constructing the n-th degree Newton polynomial 
\(P \((x)\)\  = \ 
        a\_0 + \(a\_1\) \((x - x\_0)\) + 
          \(a\_2\) \((x - x\_0)\)\ \((x - x\_1)\) + 
          \(a\_3\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  + ⋯ + 
          \(a\_n\) \((x - x\_0)\)\ \((x - x\_1)\)\ \((x - x\_2)\)\ 
            ⋯ \((x - x\_\(n - 1\))\)\)
and using  f'((Subscript[x, 0]))≈P'((Subscript[x, 0]))  as the final answer.  The method must be used at  (Subscript[x, 0]).  The points can be rearranged  {(Subscript[x, k]), (Subscript[x, 0]), (Subscript[x, 1]), ⋯, (Subscript[x, k - 1]), 
        ⋯, (Subscript[x, k + 1]), ⋯, (Subscript[x, n])}  to compute  f'((Subscript[x, k]))≈P'((Subscript[x, k])).
Section  6.2   Numerical Differentiation Formulas   Page  342

Execute this cell to activate Differentiation Based on n+1 Nodes.
DerAtX1[XY0_]
XY0_ is the list of points.

Off[General::"spell1"]; Clear[DerAtX1, XY0]; \(DerAtX1[XY0_] := 
	
      Module[{j, k, n, X, Y, XY}, 
		XY = XY0; 
		n = Length[XY] - 1; 

        		X = \(Transpose[XY]\)[[1]]; 
        
		Y = 
          \(Transpose[XY]\)[[2]]; 
		
        a = Table[0, {n + 1}]; 
		
        For[k = 0, k ≤ n, \(k++\), 
			
          a[[k + 1]] = 
            Y[[k + 1]]]; 
		
        For[j = 1, j ≤ n, \(j++\), 
			
          For[k = n, j ≤ k, \(k--\), 
            a[[k + 1]] = 
              \(a[[k + 1]] - 
                  a[[k]]\)\/\(X
                    [[k + 1]] - 
                  X[[k - j + 1]]\)]]; 

        		For[k = 0, k ≤ n, \(k++\), 
          a[[k + 1]] = 
            Together[a[[k + 1]]]]; 
		
        z = X[[1]]; 
		
        Df = a[[2]]; 
		Pr = 1; 
		
        For[k = 2, k ≤ n, \(k++\), 
			
          Pr = Pr\ \((z - X[[k]])\); 

          			Df = 
            Df + Pr\ a[[k + 1]]]; 
		
        Return[Df]]; 
\(On[General::"spell1"]; \)\)

Chapter 6, Problem 9, Page 344.
Approximate the derivative at each of the four points given in the table:

XY
       = {{0.0, 0.989992}, {0.1, 0.999135}, {0.2, 0.998295}, {0.3, 
          0.987480}};

At  x = 0.0

xy  =  XY;
x0  = xy[[1, 1]];
y0  = xy[[1, 2]];
d0  =  DerAtX1[xy];
dots = ListPlot[xy, PlotStyle -> {PointSize[0.025]}, 
          DisplayFunction -> Identity];
\(line = Graphics[{RGBColor[1, 0, 0], 
			
            Line[{{x0 - 1\/20, y0 - d0\/20}, {x0 + 1\/20, y0 + d0\/20}}]}]\);
Show[dots, line, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y    "}];
Print["
f'[",xy[[1, 1]],"] = ",d0];

Reorder the list.
At  x = 0.1

xy = {XY[[2]], 
          XY[[1]], 
          XY[[3]], 
          XY[[4]]};
x0  = xy[[1, 1]];
y0  = xy[[1, 2]];
d0  =  DerAtX1[xy];
dots = ListPlot[xy, PlotStyle -> {PointSize[0.025]}, 
          DisplayFunction -> Identity];
\(line = Graphics[{RGBColor[1, 0, 0], 
			
            Line[{{x0 - 1\/20, y0 - d0\/20}, {x0 + 1\/20, y0 + d0\/20}}]}]\);
Show[dots, line, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y    "}];
Print["
f'[",xy[[1, 1]],"] = ",d0];

Reorder the list.
At  x = 0.2

xy = {XY[[3]], 
          XY[[2]], 
          XY[[1]], 
          XY[[4]]};
x0  = xy[[1, 1]];
y0  = xy[[1, 2]];
d0  =  DerAtX1[xy];
dots = ListPlot[xy, PlotStyle -> {PointSize[0.025]}, 
          DisplayFunction -> Identity];
\(line = Graphics[{RGBColor[1, 0, 0], 
			
            Line[{{x0 - 1\/20, y0 - d0\/20}, {x0 + 1\/20, y0 + d0\/20}}]}]\);
Show[dots, line, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"x", "y    "}];
Print["
f'[",xy[[1, 1]],"] = ",d0];

Reorder the list.
At x = 0.3

xy = {XY[[4]], 
            XY[[3]], 
            XY[[2]], 
            XY[[1]]};x0  = xy[[1, 1]];
y0  = xy[[1, 2]];
d0  =  DerAtX1[xy];
dots = 
          ListPlot[xy, PlotStyle -> {PointSize[0.025]}, 
            DisplayFunction -> Identity];
\(line = 
          Graphics[{RGBColor[1, 0, 0], 
			
              Line[{{x0 - 1\/20, y0 - d0\/20}, {x0 + 1\/20, y0 + d0\/20}}]}]
            \);
Show[dots, line, DisplayFunction -> $DisplayFunction, 
          AxesLabel -> {"x", "y    "}];
Print["
f'[",xy[[1, 1]],"] = ",d0];