Module for Derivation of Numerical Differentiation Formulae

Numerical Methods using Mathematica  
Complementary software supplements for   
Numerical Methods textbooks  
John H. Mathews and Kurtis Fink, 2004

Computer Derivations of Numerical Differentiation Formulae (Classroom Notes)
John Mathews
Int. J. of Math. Education in Sci. and Tech., V 34, No 2 (March-April 2003), pp.280-287.

The traditional "pencil and paper" derivations for numerical differentiation formulas for  f'[(Subscript[x, 0])] and  f''[(Subscript[x, 0])] are done independently as if there was no connection among the derivations.  This new approach gives a parallel development of the formulas.  It requires the solution of a "linear system" that includes symbolic quantities as coefficients and constants.  The power of a computer algebra system such as Mathematica is used to elegantly solve the linear system for  f'[(Subscript[x, 0])] and  f''[(Subscript[x, 0])].  The extension to proofs of higher order numerical differentiation formulas for the central, forward or backward differences is easy to accomplish.  The notation convention using square brackets "[" and "]" instead of parenthesis is used to make all the formulas look compatible with the Mathematica input and output that is used.

The Three Point Central Difference Formulas

	Using three points f[(Subscript[x, 0]) - h], f[(Subscript[x, 0])], and f[(Subscript[x, 0]) + h]we give a parallel development of the numerical differentiation formulas for  f'[(Subscript[x, 0])] and  f''[(Subscript[x, 0])].  Start with the Taylor series for f[(Subscript[x, 0]) + h] expanded in powers of  h.  Terms must be included so that the remainder term is included.  Using degree n = 4 will suffice.

n = 4;   
s[h_] = Series[f[(Subscript[x, 0]) + h], {h, 0, n}];
\(\(Print[\*""f[\!\(x\_0\)+h] = "", s[h]\ ];\)\)

The Mathematica notation for the remainder term is  O[h]^(n + 1)=O[h]^5.  The term O[h]^5 is removed from the series with the command Normal, and a Taylor polynomial of degree n = 4 is formed.  However, we must remember that the accuracy for the Taylor polynomial is  O[h]^5.

s[h_] = 
        Normal[Series[
            f[(Subscript[x, 0]) + h], {h, 0, n}]];  
\(\(Print[\*""f[\!\(x\_0\)+h] ≈ "", 
        s[h]\ ];\)\)

Construct two "equations" by setting the series equal to the function at x = (Subscript[x, 0]) ± h.

eqns = {s[- h] == f[(Subscript[x, 0]) -  h], 
          s[ h] == f[(Subscript[x, 0]) +  h]};   
Print[ TableForm[eqns] ];

We have two equations in the two unknown f'[(Subscript[x, 0])] and  f''[(Subscript[x, 0])], and all the other quantities h,  f[(Subscript[x, 0]) - h], f[(Subscript[x, 0])], f[(Subscript[x, 0]) + h], Derivative[3][f][(Subscript[x, 0])] and  Derivative[4][f][(Subscript[x, 0])] are considered as constants.  Now use the Mathematica command "Solve" to solve for  f'[(Subscript[x, 0])] and  f''[(Subscript[x, 0])].

solset = 
        Solve[eqns, {f'[(Subscript[x, 0])], f''[
              (Subscript[x, 0])]}];   
Print["Solve the equations"]; 
Print[ TableForm[eqns] ];
Print["Get"]; 
Print[ 
        TableForm[
          (Subscript[solset, [[1]]])] ];

The result is not too easy to read, so we can use the commands ExpandAll and Collect to manipulate the formula.  In the print statements we add  "+ ... "  to remind us that we are using truncated infinite series.  The following commands will create a nicer looking printout for the above solution.

sol1 = 
        Collect[ExpandAll[
            (Subscript[solset, [[1, 1, 2]]])], h];  
(Subscript[sol1, [[1]]]) = 
        Together[
          (Subscript[sol1, [[1]]])];  
sol2 = 
        Collect[ExpandAll[
            (Subscript[solset, [[1, 2, 2]]])], h];  
\(\(Print[\*"" f'[\!\(x\_0\)] = "", 
        sol1, " + ..."];\)\ \ \ \)
\(\(Print[\*""f''[\!\(x\_0\)] = "", 
        sol2, " + ..."];\)\ \ \ \)

Thus, we have derived the numerical differentiation formula for   f'[(Subscript[x, 0])]  and   f''[(Subscript[x, 0])] and the first term in the series expansion for the remainder which involves, Derivative[3][f][(Subscript[x, 0])] or  Derivative[4][f][(Subscript[x, 0])], respectively.  Since the "numerical differentiation formulae" are "truncated" infinite series, we know that the ellipses "..." means that there are infinitely many more term which are not shown.

When we do not include the  "...", we must evaluate the lowest order derivative in the series for the remainder at the value x = c instead of x = (Subscript[x, 0]), then we can "chop off" the infinite series at the term involving  Derivative[3][f][c] or  Derivative[4][f][c], respectively.  The Mathematica command  ReplaceAll is used to accomplish this task.

sol11=ReplaceAll[Take[sol1, 2],Derivative[k_][f][(Subscript[x, 0])]->Derivative[k][f][c]];  
sol22=ReplaceAll[Take[sol2, 2],Derivative[k_][f][(Subscript[x, 0])]->Derivative[k][f][c]];  
\(Print[\*"" f'[\!\(x\_0\)] = "", sol11];\)   
\(Print[\*""f''[\!\(x\_0\)] = "", sol22];\)

We have obtained the desired numerical differentiation formulas and their remainder terms.  We should mention that new in version 4.1 of  Mathematica is symbolic recognition of order of approximation for divided difference formulas.  We can add the term O[h]^2 to the numerical differentiation formula term, and Mathematica will conclude that it is the derivative plus O[h]^2.   The following commands illustrate this feature of Mathematica.

\(\(\(Print[\*""((-f[-h + (Subscript[x, 0])] + f[h + (Subscript[x, 0])])/(2  h))+ O[h]^2 = "", \(\(-f[\(-h\) + x\_0]\) + f[h + x\_0]\)\/\(2\ h\) + \ O[h]\^2];\)\(\ \)\)\)

\(\(Print[\*""((-2 f[(Subscript[x, 0])] + f[-h + (Subscript[x, 0])] + f[h + (Subscript[x, 0])])/(h^2))+ O[h]^2 = "", \(\(-2\)\ f[x\_0] + f[\(-h\) + x\_0] + f[h + x\_0]\)\/h\^2 + \ O[h]\^2];\)\)

Therefore, we have established the numerical differentiation formulas   

	  f'[(Subscript[x, 0])]= ((f[(Subscript[x, 0]) + h]  -  f[(Subscript[x, 0]) - h])/(2 h)) - (1/6) Derivative[3][f][c](h^2)  
and
	f''[(Subscript[x, 0])]= ((f[(Subscript[x, 0]) - h]  -  2 f[(Subscript[x, 0])]  +  f[(Subscript[x, 0]) + h])/(h^2)) - (1/12) Derivative[4][f][c] (h^2)

And the corresponding formulas with the big "O" notation  O[h]^2 are   

	  f'[(Subscript[x, 0])]= ((f[(Subscript[x, 0]) + h]  -  f[(Subscript[x, 0]) - h])/(2 h)) + O[h]^2  
and
	f''[(Subscript[x, 0])]=((f[(Subscript[x, 0]) - h]  -  2 f[(Subscript[x, 0])]  +  f[(Subscript[x, 0]) + h])/(h^2)) + O[h]^2

Comparison with the Traditional Derivations

	The above derivation differs only slightly from the traditional derivations which also start with the two equations  

	f[(Subscript[x, 0]) - h] = f[(Subscript[x, 0])] - h f'[(Subscript[x, 0])] + (1/2) (h^2) f''[(Subscript[x, 0])] - (1/6) (h^3) Derivative[3][f][(Subscript[x, 0])]+(1/24) (h^4) Derivative[4][f][(Subscript[x, 0])]+ ...  
(1) 
	f[h + (Subscript[x, 0])] = f[(Subscript[x, 0])] + h f'[(Subscript[x, 0])] + (1/2) (h^2) f''[(Subscript[x, 0])] + (1/6) (h^3) Derivative[3][f][(Subscript[x, 0])]+(1/24) (h^4) Derivative[4][f][(Subscript[x, 0])]+ ...

When deriving the numerical differentiation formula for f'[(Subscript[x, 0])] the equations (1) are subtracted and the terms involving  f''[(Subscript[x, 0])] cancel and manipulations are used to solve for  f'[(Subscript[x, 0])] and its truncation error term.

The traditional derivation of the numerical differentiation formula for f''[(Subscript[x, 0])] starts with the same equations (1).  But this time the equations are added and the terms involving  f'[(Subscript[x, 0])] cancel and manipulations are used to solve for  f''[(Subscript[x, 0])] and its truncation error term.

There is no "big deal" made about the fact that the starting place is the same.  We shall see for the higher order formulas that using the same starting place will be the key to successful computer derivations of numerical differentiation formulas.

The Five Point Central Difference Formulas

	Using five points f[(Subscript[x, 0]) - 2  h], f[(Subscript[x, 0]) - h], f[(Subscript[x, 0])],f[(Subscript[x, 0]) + h], and f[(Subscript[x, 0]) + 2  h]we can give a parallel development of the numerical differentiation formulas for  f'[(Subscript[x, 0])],  f''[(Subscript[x, 0])], Derivative[3][f][(Subscript[x, 0])] and Derivative[4][f][(Subscript[x, 0])].  Start with the Taylor series for f[(Subscript[x, 0]) + h] expanded in powers of  h.  Terms must be included so that the remainder term is included.  Using degree n = 6 will suffice.  Since the remainder term  O[h]^(n + 1)=O[h]^7, we will use Mathematica's Normal command to remove it from the series, and form the Taylor polynomial of degree n = 6.

k = 2;  n = 2  k + 2;   
s[h_] = 
        Normal[Series[
            f[(Subscript[x, 0]) + h], {h, 0, n}]];  
\(\(Print[\*""f[\!\(x\_0\)+h] ≈ "", 
        s[h]\ ];\)\)

Use five points and set up the "equations" by setting  s[(Subscript[x, 0]) + k h] == f[(Subscript[x, 0]) + k h] for  k = -2,-1,0,1,2, which is automated by using Mathematica's Table command.  Although the value   k = 0, will produce  f[(Subscript[x, 0])] == f[(Subscript[x, 0])] which is the boolean value True, it will do no harm when solving the set of equations and we will see that the generations of higher order formulas which require more complicated sets of equations is more easily automated with the use of a table.

eqns = 
        Table[s[j h] == f[(Subscript[x, 0]) +  j h], {j, -k, 
            k}];  
Print[ TableForm[eqns] ];

We can consider that these are four equations in the four unknowns f'[(Subscript[x, 0])], f''[(Subscript[x, 0])],  Derivative[3][f][(Subscript[x, 0])] and Derivative[4][f][(Subscript[x, 0])], and all the other quantities h, f[(Subscript[x, 0]) - 2  h], f[(Subscript[x, 0]) - h], f[(Subscript[x, 0])], f[(Subscript[x, 0]) + h], f[(Subscript[x, 0]) - 2  h], Derivative[5][f][(Subscript[x, 0])] and  Derivative[6][f][(Subscript[x, 0])] are constants.  Now use Mathematica to solve for f'[(Subscript[x, 0])], f''[(Subscript[x, 0])],  Derivative[3][f][(Subscript[x, 0])] and Derivative[4][f][(Subscript[x, 0])].  Since this requires typing of four derivatives to be solved, we will automate this process too, by using a table to construct the "variables."

vars = Table[Derivative[k][f][(Subscript[x, 0])],{k, 1, 2  k}];  
Print["The variables in the above \
equations are"];
Print[ vars ];

Amazingly, the algebra involved in solving the four equations will result in the cancellation of Derivative[6][f][(Subscript[x, 0])]in the odd order derivatives and the term  Derivative[5][f][(Subscript[x, 0])] will cancel in the even order derivatives.  This can be accomplished using Mathematica's Solve procedure.

solset = Solve[eqns, vars];   
Print["Solve the equations"]; 
Print[ TableForm[eqns] ];
Print["Get"]; 
Print[ 
        TableForm[
          (Subscript[solset, [[1]]])] ];

The output is not too clear to read, so we shall use the formula manipulation commands ExpandAll and Collect to group the numerical differentiation portion, and then use the ReplaceAll command to change Derivative[5][f][(Subscript[x, 0])] and  Derivative[6][f][(Subscript[x, 0])] to Derivative[5][f][c] and  Derivative[6][f][c], then a final use of Together will clean up the numerical differentiation formula part.  All this will be accomplished in the following three Mathematica commands.  If the reader is curious about exactly what is happening in each step, then the semi-colons can be deleted and the results of each operation will be shown.

(Subscript[solset, [[1, All, 2]]]) = 
          Collect[ExpandAll[
              (Subscript[solset, [[1, All, 2]]])], h];  
(Subscript[solset, [[1, All, 2, 
              2]]])=ReplaceAll[(Subscript[solset, [[1, All, 2, 
                  2]]]),Derivative[i_][f][(Subscript[x, 0])]->Derivative[i][f][c]];  
\(solset\_\(\([[\)\(1, All, 2, 1\)\(]]\)\) = 
        Together[
          solset\_\(\([[\)\(1, All, 2, 1\)\(]]\)\)]; \ 
      Print[\*"" f'[(Subscript[x, 0])] = "", 
        solset\_\(\([[\)\(1, 1, 2\)\(]]\)\)];\)   
\(Print[\*""f''[(Subscript[x, 0])] = "", 
          solset\_\(\([[\)\(1, 2, 2\)\(]]\)\)];\)   
Print["
(f^((3)))[(Subscript[x, 0])] = ",(Subscript[solset, [[1, 3, 
                2]]])];   
Print["
(f^((4)))[(Subscript[x, 0])] = ",(Subscript[solset, [[1, 4, 
                2]]])];

Notice that the formula for approximating f'[(Subscript[x, 0])] and f''[(Subscript[x, 0])] have truncation error terms involving  (h^4) so they are numerical differentiation formulas of order O[h]^4.   But the formula for approximating Derivative[3][f][(Subscript[x, 0])] and Derivative[4][f][(Subscript[x, 0])] have truncation error terms involving  (h^2) so they are numerical differentiation formulas of order O[h]^2.  This is one of the many surprises in the theory of numerical analysis.

Again, we can uses Mathematica's 4.1 ability to symbolic recognize the order of approximation for divided difference formulas.  If we add the term  O[h]^4 or O[h]^2 to the numerical differentiation formula then we will be surprised to learn that Mathematica will tell us that it is the derivative plus  O[h]^4 or O[h]^2 .   The following commands illustrate this point.

\(\(\(Print[\*""((f[-2 h + (Subscript[x, 0])] - 8  f[-h + (Subscript[x, 0])] + 8  f[h + (Subscript[x, 0])] - f[2  h + (Subscript[x, 0])])/(12  h))+ O[h]^4 = "", \(f[\(-2\)\ h + x\_0] - 8\ f[\(-h\) + x\_0] + 8\ f[h + x\_0] - f[2\ h + x\_0]\)\/\(12\ h\) + \ O[h]\^4];\)\(\ \)\)\)

\(\(\(Print[\*""((-30 f[(Subscript[x, 0])] - f[-2 h + (Subscript[x, 0])] + 16  f[-h + (Subscript[x, 0])] + 16  f[h + (Subscript[x, 0])] - f[2  h + (Subscript[x, 0])])/(12  (h^2)))+ O[h]^4 = "", \(\(-30\)\ f[x\_0] - f[\(-2\)\ h + x\_0] + 16\ f[\(-h\) + x\_0] + 16\ f[h + x\_0] - f[2\ h + x\_0]\)\/\(12\ h\^2\) + \ 
          O[h]\^4];\)\(\ \)\)\)

\(\(\(Print[\*""((-f[-2 h + (Subscript[x, 0])] + 2  f[-h + (Subscript[x, 0])] - 2  f[h + (Subscript[x, 0])] + f[2  h + (Subscript[x, 0])])/(2  (h^3)))+ O[h]^2 = "", \(\(-f[\(-2\)\ h + x\_0]\) + 2\ f[\(-h\) + x\_0] - 2\ f[h + x\_0] + f[2\ h + x\_0]\)\/\(2\ h\^3\) + \ O[h]\^2];\)\(\ \)\)\)

\(\(\(Print[\*""((6  f[(Subscript[x, 0])] + f[-2 h + (Subscript[x, 0])] - 4  f[-h + (Subscript[x, 0])] - 4  f[h + (Subscript[x, 0])] + f[2  h + (Subscript[x, 0])])/(h^4)) +O[h]^2 = "", \(6\ f[x\_0] + f[\(-2\) h + x\_0] - 4  f[\(-h\) + x\_0] - 4  f[h + x\_0] + f[2  h + x\_0]\)\/h\^4\  + O[h]\^2];\)\(\ \)\)\)

Therefore, we have established the numerical differentiation formulas   

	 f'[(Subscript[x, 0])]= \(\(f[x\_0 - 2\ h]\  - \ 8\ f[x\_0 - h]\  + \ 8\ f[x\_0 + h]\  - \ f[x\_0 + 2\ h]\)\/\(12\ h\)\) +(1/30) (h^4) Derivative[5][f][c]    


	f''[(Subscript[x, 0])]= \(\(\(-f[x\_0 - 2  h]\)\  + \ 16\ f[x\_0 - h]\  - \ 30\ f[x\_0]\  + \ 16\ f[x\_0 + h]\  - \ f[x\_0 + 2  h]\)\/\(12\ h\^2\)\) +(1/90) (h^4) Derivative[6][f][c]  


	Derivative[3][f][(Subscript[x, 0])]= \(\(\(-f[x\_0 - 2  h]\)\  + \ 2\ f[x\_0 - h]\  - \ 2\ f[x\_0 + h]\  + \ f[x\_0 + 2  h]\)\/\(2\ h\^3\)\) -(1/4) (h^2) Derivative[5][f][c]   


	Derivative[4][f][(Subscript[x, 0])]= \(\(f[x\_0 - 2  h]\  - \ 4\ f[x\_0 - h]\  + \ 6\ f[x\_0]\  - \ 4\ f[x\_0 + h]\  + \ f[x\_0 + 2  h]\)\/h\^4\) -(1/6) (h^2) Derivative[6][f][c]

These formulas can be written with the big "O" notation O[h]^4 and O[h]^2 if desired.

	 f'[(Subscript[x, 0])]= \(\(f[x\_0 - 2\ h]\  - \ 8\ f[x\_0 - h]\  + \ 8\ f[x\_0 + h]\  - \ f[x\_0 + 2\ h]\)\/\(12\ h\)\) + O[h]^4    


	f''[(Subscript[x, 0])]= \(\(\(-f[x\_0 - 2  h]\)\  + \ 16\ f[x\_0 - h]\  - \ 30\ f[x\_0]\  + \ 16\ f[x\_0 + h]\  - \ f[x\_0 + 2  h]\)\/\(12\ h\^2\)\) + O[h]^4   


	Derivative[3][f][(Subscript[x, 0])]= \(\(\(-f[x\_0 - 2  h]\)\  + \ 2\ f[x\_0 - h]\  - \ 2\ f[x\_0 + h]\  + \ f[x\_0 + 2  h]\)\/\(2\ h\^3\)\) + O[h]^2   


	Derivative[4][f][(Subscript[x, 0])]= \(\(f[x\_0 - 2  h]\  - \ 4\ f[x\_0 - h]\  + \ 6\ f[x\_0]\  - \ 4\ f[x\_0 + h]\  + \ f[x\_0 + 2  h]\)\/h\^4\) + O[h]^2

Appendix. Subroutines for Generating Numerical Differentiation Formulae

The central difference formulae use an odd number of points  m = 2  k + 1,  and an even number of  n = 2  k equations.  Since we want to include the remainder terms, we need to use series expansions of order  n = 2  k + 2.  The remainder term in these formulas all involve even powers of  h  and even and odd derivatives depending on the situation.  Also, the subroutine requires a replacement of the point where the remainder term is evaluated to be x = c instead of  x = (Subscript[x, 0]).

CentralDiff[k0_]:=Module[{k = k0},
Clear[h, i, j, n, f, s, x];  
s[h_] = 
                  Normal[Series[f[(Subscript[x, 0]) + h], {h, 0, 2  k + 2}]];  
eqns = 
                  Table[s[j h] == f[(Subscript[x, 0]) +  j h], {j, -k, k}];    
vars = Table[Derivative[k][f][(Subscript[x, 0])],{k, 1, 2  k}];  
pts  =  
                  Table[{(Subscript[x, 0]), f[(Subscript[x, 0]) +  j h]}, {j, -k, k}];  
solset = Solve[eqns, vars];  
(Subscript[solset, [[1, All, 2]]]) = Collect[ExpandAll[
                      (Subscript[solset, [[1, All, 2]]])], h];  
(Subscript[solset, [[1, All, 2, 
                      2]]])=ReplaceAll[(Subscript[solset, [[1, All, 2, 
                          2]]]),Derivative[i_][f][(Subscript[x, 0])]->Derivative[i][f][c]];  
(Subscript[solset, [[1, All, 2, 1]]]) = Together[
                    (Subscript[solset, [[1, All, 2, 1]]])];  
Print["Central difference formulas \
for numerical differentiation"];
Print["Using the ", 
                  2  k + 1, " points"];
Print[ pts, "\n" ];
Print[
                  TableForm[
                    (Subscript[solset, [[1]]])]]; ];

The Five Point Central Difference Formulae

	These are to difficult to do by hand, so we use the computer.

CentralDiff[2]

Higher order formulas are easy to obtain with Mathematica.

Exploration

CentralDiff[3]

CentralDiff[4]

Subroutines for Generating the Forward and Backward Difference Formulae

Once we have the idea for the central difference formula, all we need to do is adjust the set of to be forward step sizes or backward step sizes and we can generate the forward and backward differentiation formulas.

ForwardDiff[k0_]:=Module[{k = k0},
Clear[h, i, j, n, f, s, x];  
s[h_] = 
                  Normal[Series[f[(Subscript[x, 0]) + h], {h, 0, k + 1}]];  
eqns = 
                  Table[s[j h] == f[(Subscript[x, 0]) +  j h], {j, 0, k}];    
vars = Table[Derivative[k][f][(Subscript[x, 0])],{k, 1, k}];  
pts  =  
                  Table[{(Subscript[x, 0]), f[(Subscript[x, 0]) +  j h]}, {j, 0, k}];  
solset = Solve[eqns, vars];  
(Subscript[solset, [[1, All, 2]]]) = 
                  Collect[
                    ExpandAll[
                      (Subscript[solset, [[1, All, 2]]])], h];  
(Subscript[solset, [[1, All, 2, 
                      2]]])=ReplaceAll[(Subscript[solset, [[1, All, 2, 
                          2]]]),Derivative[i_][f][(Subscript[x, 0])]->Derivative[i][f][c]];  
(Subscript[solset, [[1, All, 2, 1]]]) = Together[
                    (Subscript[solset, [[1, All, 2, 1]]])];  
Print["Forward difference formulas \
for numerical differentiation"];
Print["Using the ", 
                  k + 1, " points"];
Print[ pts, "\n" ];
Print[
                  TableForm[
                    (Subscript[solset, [[1]]])]]; ];

BackwardDiff[k0_]:=Module[{k = k0},
Clear[h, j, n, f, s, x];  
s[h_] = 
                  Normal[Series[f[(Subscript[x, 0]) + h], {h, 0, k + 1}]];  
eqns = 
                  Table[s[-j h] == f[(Subscript[x, 0]) -  j h], {j, 0, k}];    
vars = Table[Derivative[k][f][(Subscript[x, 0])],{k, 1, k}];  
pts  =  
                  Table[{(Subscript[x, 0]), f[(Subscript[x, 0]) +  j h]}, {j, 0, k}];  
solset = Solve[eqns, vars];  
(Subscript[solset, [[1, All, 2]]]) = 
                  Collect[
                    ExpandAll[
                      (Subscript[solset, [[1, All, 2]]])], h];  
(Subscript[solset, [[1, All, 2, 
                      2]]])=ReplaceAll[(Subscript[solset, [[1, All, 2, 
                          2]]]),Derivative[i_][f][(Subscript[x, 0])]->Derivative[i][f][c]];  
(Subscript[solset, [[1, All, 2, 1]]]) = Together[
                    (Subscript[solset, [[1, All, 2, 1]]])];  
Print["Backward difference \
formulas for numerical differentiation"];
Print["Using the ", 
                  k + 1, " points"];
Print[ pts, "\n" ];
Print[
                  TableForm[
                    (Subscript[solset, [[1]]])]]; ];

Explorations with the Higher Order Formulas

BackwardDiff[3]

ForwardDiff[3]

Higher order formulas are easy to obtain with Mathematica.

Exploration

BackwardDiff[1]

BackwardDiff[2]

BackwardDiff[3]

BackwardDiff[4]

BackwardDiff[5]

BackwardDiff[6]

BackwardDiff[7]

Reference

John Mathews, Computer Derivations of Numerical Differentiation Formulae (Classroom Notes),  International Journal of Mathematics Education in Science and Technology, Volume 34 No 2 (March-April 2003), pp.280-287.

Research Experience for Undergraduates

Numerical Differentiation  Numerical Differentiation  Internet hyperlinks to web sites and a bibliography of articles.