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

Algorithm 6.1 (Differentiation Using Limits).  
To approximate f'(x)  numerically by generating the sequence 

			f'(x)  ~  dk  =  ( f(x+h/2-k)Ê-Êf(x-h/2-k) )/( h/2-k+1 )   
			
for  k = 0,1,...,n   until   |dn+1 - dn| <= |dn - dn-1|   or   |dn - dn-1| < Tol,   
which is an attempt to find the best  approximation   f'(x)  ~  dn.
Section	6.1	Approximating the Derivative   Page  326

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 = 10^-7 //N;
  big = 10^6 //N;
  max = 32;
  min = 2;
  Clear[d,e];
  d = {};
  e = {1};
  h = h0 //N;
  dydx = 0.5 (f[x0+h] - f[x0-h])/h //N;
  d = Append[d,dydx];
  n = 2;
  For[j=2,j<=3,j++,
    h = h/2;
    dydx = 0.5 (f[x0+h] - f[x0-h])/h //N;
    Print[N[dydx,14],"\t     h = ",h];
    d = Append[d,dydx];
    n = Length[d];
    Err = Abs[d[[n]] - d[[n-1]]];
    e = Append[e,Err];
    RelErr=2Err/(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 = 0.5 (f[x0+h] - f[x0-h])/h //N;
    Print[N[dydx,14],"\t     h = ",h];
    d = Append[d,dydx];
    n = Length[d];
    Err = Abs[d[[n]] - d[[n-1]]];
    e = Append[e,Err];
    RelErr=2Err/(Abs[d[[n]]]+Abs[d[[n-1]]]+small) ];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to f'[x] 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]

Compare with the analytic method of computing derivatives.

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) )/( 4kÊ-Ê1 )    for 1 <= k <=  j.
Section	6.1	Approximating the Derivative   Page  327

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 = 10^-7 //N;
  big = 10^6 //N;
  max = 32;
  min = 2;
  Clear[d,e];
  d = Table[0,{11},{11}];
  h = h0 //N;
  j = 1;
  Error = 1.0;
  RelErr = 1.0;
  d[[1,1]] = 0.5 (f[x0+h] - f[x0-h])/h //N;
  Print[N[d[[1,1]],14],"\t     h = ",h];
  While[RelErr>Tol && Error>Epsilon && j<10,
    h = h/2;
    d[[j+1,1]] = 0.5 (f[x0+h] - f[x0-h])/h //N;
    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],"\t     h = ",h];
    Error = Abs[d[[j+1,j+1]]-d[[j,j]]] //N;
    RelErr = 2Error/(Abs[d[[j+1,j+1]]] +
             Abs[d[[j,j]]] + small) //N;
    j = j+1;
    n = j; ];
  Print[" "];
  Print["The function is f[x] = ",f[x]];
  Print["A numerical approximation to f'[x] 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]

Compare with the analytic method of computing derivatives.

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)  =  a0 + a1(x - x0) + a2(x - x0)(x - x1) + a3(x - x0)(x - x1)(x - x2)
			      +... + an(x - x0)(x - x1)...(x - xn-1),
			  
and using  f'(x0) ~ P'(x0)  as the final answer.  
The method must be used at x0.  
The points can be rearranged  {xk,x0,...,xk-1,xk+1,...,xn}  to compute         
f'(xk) ~ P'(xk).
Section	6.2	Numerical Differentiation Formulas   Page  342

Off[General::spell1];
Clear[DerAtX1,XY0];
DerAtX1[XY0_] := 
Block[{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

DerAtX1[XY]

Reorder the list.

At x = 0.1

xy = {XY[[2]],XY[[1]],XY[[3]],XY[[4]]}

DerAtX1[xy]

Reorder the list.

At x = 0.2

xy = {XY[[3]],XY[[2]],XY[[1]],XY[[4]]}

DerAtX1[xy]

Reorder the list.

At x = 0.3

xy = {XY[[4]],XY[[3]],XY[[2]],XY[[1]]}

DerAtX1[xy]