Module

for

Derivation of Numerical Differentiation Formulae

     

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  [Graphics:Images/NumericalDiffFormulaeMod_gr_1.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_2.gif] 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  [Graphics:Images/NumericalDiffFormulaeMod_gr_3.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_4.gif].  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 [Graphics:Images/NumericalDiffFormulaeMod_gr_5.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_6.gif], and [Graphics:Images/NumericalDiffFormulaeMod_gr_7.gif]we give a parallel development of the numerical differentiation formulas for  [Graphics:Images/NumericalDiffFormulaeMod_gr_8.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_9.gif].  Start with the Taylor series for [Graphics:Images/NumericalDiffFormulaeMod_gr_10.gif] expanded in powers of  [Graphics:Images/NumericalDiffFormulaeMod_gr_11.gif].  Terms must be included so that the remainder term is included.  Using degree [Graphics:Images/NumericalDiffFormulaeMod_gr_12.gif] will suffice.    

[Graphics:Images/NumericalDiffFormulaeMod_gr_13.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_14.gif]

The Mathematica notation for the remainder term is  [Graphics:Images/NumericalDiffFormulaeMod_gr_15.gif].  The term [Graphics:Images/NumericalDiffFormulaeMod_gr_16.gif] is removed from the series with the command Normal, and a Taylor polynomial of degree [Graphics:Images/NumericalDiffFormulaeMod_gr_17.gif] is formed.  However, we must remember that the accuracy for the Taylor polynomial is  [Graphics:Images/NumericalDiffFormulaeMod_gr_18.gif].   

[Graphics:Images/NumericalDiffFormulaeMod_gr_19.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_20.gif]

Construct two "equations" by setting the series equal to the function at [Graphics:Images/NumericalDiffFormulaeMod_gr_21.gif].    

[Graphics:Images/NumericalDiffFormulaeMod_gr_22.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_23.gif]

 

We have two equations in the two unknown [Graphics:Images/NumericalDiffFormulaeMod_gr_24.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_25.gif], and all the other quantities [Graphics:Images/NumericalDiffFormulaeMod_gr_26.gif],  [Graphics:Images/NumericalDiffFormulaeMod_gr_27.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_28.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_29.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_30.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_31.gif] are considered as constants.  Now use the Mathematica command "Solve" to solve for  [Graphics:Images/NumericalDiffFormulaeMod_gr_32.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_33.gif].    

[Graphics:Images/NumericalDiffFormulaeMod_gr_34.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_35.gif]

The result is not too easy to read, so we can use the commands [Graphics:Images/NumericalDiffFormulaeMod_gr_36.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_37.gif] to manipulate the formula.  In the print statements we add  "[Graphics:Images/NumericalDiffFormulaeMod_gr_38.gif]"  to remind us that we are using truncated infinite series.  The following commands will create a nicer looking printout for the above solution.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_39.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_40.gif]

 

Thus, we have derived the numerical differentiation formula for   [Graphics:Images/NumericalDiffFormulaeMod_gr_41.gif]  and   [Graphics:Images/NumericalDiffFormulaeMod_gr_42.gif] and the first term in the series expansion for the remainder which involves, [Graphics:Images/NumericalDiffFormulaeMod_gr_43.gif] or  [Graphics:Images/NumericalDiffFormulaeMod_gr_44.gif], respectively.  Since the "numerical differentiation formulae" are "truncated" infinite series, we know that the ellipses "[Graphics:Images/NumericalDiffFormulaeMod_gr_45.gif]" means that there are infinitely many more term which are not shown.   

    When we do not include the  "[Graphics:Images/NumericalDiffFormulaeMod_gr_46.gif]", we must evaluate the lowest order derivative in the series for the remainder at the value [Graphics:Images/NumericalDiffFormulaeMod_gr_47.gif] instead of [Graphics:Images/NumericalDiffFormulaeMod_gr_48.gif], then we can "chop off" the infinite series at the term involving  [Graphics:Images/NumericalDiffFormulaeMod_gr_49.gif] or  [Graphics:Images/NumericalDiffFormulaeMod_gr_50.gif], respectively.  The Mathematica command  ReplaceAll is used to accomplish this task.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_51.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_52.gif]

 

    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 [Graphics:Images/NumericalDiffFormulaeMod_gr_53.gif] to the numerical differentiation formula term, and Mathematica will conclude that it is the derivative plus [Graphics:Images/NumericalDiffFormulaeMod_gr_54.gif].   The following commands illustrate this feature of Mathematica.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_55.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_56.gif]

   

 

[Graphics:Images/NumericalDiffFormulaeMod_gr_57.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_58.gif]

Therefore, we have established the numerical differentiation formulas   

        [Graphics:Images/NumericalDiffFormulaeMod_gr_59.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_60.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_61.gif]  
and
        [Graphics:Images/NumericalDiffFormulaeMod_gr_62.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_63.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_64.gif]

And the corresponding formulas with the big "O" notation  [Graphics:Images/NumericalDiffFormulaeMod_gr_65.gif] are   

        [Graphics:Images/NumericalDiffFormulaeMod_gr_66.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_67.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_68.gif]  
and
        [Graphics:Images/NumericalDiffFormulaeMod_gr_69.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_70.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_71.gif]

Comparison with the Traditional Derivations

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

        [Graphics:Images/NumericalDiffFormulaeMod_gr_72.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_73.gif] [Graphics:Images/NumericalDiffFormulaeMod_gr_74.gif]  
(1)
        [Graphics:Images/NumericalDiffFormulaeMod_gr_75.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_76.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_77.gif]  

    When deriving the numerical differentiation formula for [Graphics:Images/NumericalDiffFormulaeMod_gr_78.gif] the equations (1) are subtracted and the terms involving  [Graphics:Images/NumericalDiffFormulaeMod_gr_79.gif] cancel and manipulations are used to solve for  [Graphics:Images/NumericalDiffFormulaeMod_gr_80.gif] and its truncation error term.  

     The traditional derivation of the numerical differentiation formula for [Graphics:Images/NumericalDiffFormulaeMod_gr_81.gif] starts with the same equations (1).  But this time the equations are added and the terms involving  [Graphics:Images/NumericalDiffFormulaeMod_gr_82.gif] cancel and manipulations are used to solve for  [Graphics:Images/NumericalDiffFormulaeMod_gr_83.gif] 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 [Graphics:Images/NumericalDiffFormulaeMod_gr_84.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_85.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_86.gif],[Graphics:Images/NumericalDiffFormulaeMod_gr_87.gif], and [Graphics:Images/NumericalDiffFormulaeMod_gr_88.gif]we can give a parallel development of the numerical differentiation formulas for  [Graphics:Images/NumericalDiffFormulaeMod_gr_89.gif],  [Graphics:Images/NumericalDiffFormulaeMod_gr_90.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_91.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_92.gif].  Start with the Taylor series for [Graphics:Images/NumericalDiffFormulaeMod_gr_93.gif] expanded in powers of  [Graphics:Images/NumericalDiffFormulaeMod_gr_94.gif].  Terms must be included so that the remainder term is included.  Using degree [Graphics:Images/NumericalDiffFormulaeMod_gr_95.gif] will suffice.  Since the remainder term  [Graphics:Images/NumericalDiffFormulaeMod_gr_96.gif], we will use Mathematica's Normal command to remove it from the series, and form the Taylor polynomial of degree [Graphics:Images/NumericalDiffFormulaeMod_gr_97.gif].   

[Graphics:Images/NumericalDiffFormulaeMod_gr_98.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_99.gif]

Use five points and set up the "equations" by setting  [Graphics:Images/NumericalDiffFormulaeMod_gr_100.gif] for  [Graphics:Images/NumericalDiffFormulaeMod_gr_101.gif], which is automated by using Mathematica's Table command.  Although the value   [Graphics:Images/NumericalDiffFormulaeMod_gr_102.gif], will produce  [Graphics:Images/NumericalDiffFormulaeMod_gr_103.gif] 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.  

[Graphics:Images/NumericalDiffFormulaeMod_gr_104.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_105.gif]

 

We can consider that these are four equations in the four unknowns [Graphics:Images/NumericalDiffFormulaeMod_gr_106.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_107.gif],  [Graphics:Images/NumericalDiffFormulaeMod_gr_108.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_109.gif], and all the other quantities [Graphics:Images/NumericalDiffFormulaeMod_gr_110.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_111.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_112.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_113.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_114.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_115.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_116.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_117.gif] are constants.  Now use Mathematica to solve for [Graphics:Images/NumericalDiffFormulaeMod_gr_118.gif], [Graphics:Images/NumericalDiffFormulaeMod_gr_119.gif],  [Graphics:Images/NumericalDiffFormulaeMod_gr_120.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_121.gif].  Since this requires typing of four derivatives to be solved, we will automate this process too, by using a table to construct the "variables."  

[Graphics:Images/NumericalDiffFormulaeMod_gr_122.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_123.gif]
[Graphics:Images/NumericalDiffFormulaeMod_gr_124.gif]

Amazingly, the algebra involved in solving the four equations will result in the cancellation of [Graphics:Images/NumericalDiffFormulaeMod_gr_125.gif]in the odd order derivatives and the term  [Graphics:Images/NumericalDiffFormulaeMod_gr_126.gif] will cancel in the even order derivatives.  This can be accomplished using Mathematica's Solve procedure.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_127.gif]


[Graphics:Images/NumericalDiffFormulaeMod_gr_128.gif]

The output is not too clear to read, so we shall use the formula manipulation commands [Graphics:Images/NumericalDiffFormulaeMod_gr_129.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_130.gif] to group the numerical differentiation portion, and then use the ReplaceAll command to change [Graphics:Images/NumericalDiffFormulaeMod_gr_131.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_132.gif] to [Graphics:Images/NumericalDiffFormulaeMod_gr_133.gif] and  [Graphics:Images/NumericalDiffFormulaeMod_gr_134.gif], 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.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_135.gif]



[Graphics:Images/NumericalDiffFormulaeMod_gr_136.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_137.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_138.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_139.gif]

Notice that the formula for approximating [Graphics:Images/NumericalDiffFormulaeMod_gr_140.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_141.gif] have truncation error terms involving  [Graphics:Images/NumericalDiffFormulaeMod_gr_142.gif] so they are numerical differentiation formulas of order [Graphics:Images/NumericalDiffFormulaeMod_gr_143.gif].   But the formula for approximating [Graphics:Images/NumericalDiffFormulaeMod_gr_144.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_145.gif] have truncation error terms involving  [Graphics:Images/NumericalDiffFormulaeMod_gr_146.gif] so they are numerical differentiation formulas of order [Graphics:Images/NumericalDiffFormulaeMod_gr_147.gif].  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  [Graphics:Images/NumericalDiffFormulaeMod_gr_148.gif] or [Graphics:Images/NumericalDiffFormulaeMod_gr_149.gif] to the numerical differentiation formula then we will be surprised to learn that Mathematica will tell us that it is the derivative plus  [Graphics:Images/NumericalDiffFormulaeMod_gr_150.gif] or [Graphics:Images/NumericalDiffFormulaeMod_gr_151.gif] .   The following commands illustrate this point.   

[Graphics:Images/NumericalDiffFormulaeMod_gr_152.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_153.gif]

   

 

[Graphics:Images/NumericalDiffFormulaeMod_gr_154.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_155.gif]

   

 

[Graphics:Images/NumericalDiffFormulaeMod_gr_156.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_157.gif]

   

 

[Graphics:Images/NumericalDiffFormulaeMod_gr_158.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_159.gif]

Therefore, we have established the numerical differentiation formulas   

        [Graphics:Images/NumericalDiffFormulaeMod_gr_160.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_161.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_162.gif]    


        [Graphics:Images/NumericalDiffFormulaeMod_gr_163.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_164.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_165.gif]  


        [Graphics:Images/NumericalDiffFormulaeMod_gr_166.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_167.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_168.gif]   


        [Graphics:Images/NumericalDiffFormulaeMod_gr_169.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_170.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_171.gif]   

These formulas can be written with the big "O" notation [Graphics:Images/NumericalDiffFormulaeMod_gr_172.gif] and [Graphics:Images/NumericalDiffFormulaeMod_gr_173.gif] if desired.

        [Graphics:Images/NumericalDiffFormulaeMod_gr_174.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_175.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_176.gif]    


        [Graphics:Images/NumericalDiffFormulaeMod_gr_177.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_178.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_179.gif]   


        [Graphics:Images/NumericalDiffFormulaeMod_gr_180.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_181.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_182.gif]   


        [Graphics:Images/NumericalDiffFormulaeMod_gr_183.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_184.gif][Graphics:Images/NumericalDiffFormulaeMod_gr_185.gif]   

 

Appendix. Subroutines for Generating Numerical Differentiation Formulae

The central difference formulae use an odd number of points  [Graphics:Images/NumericalDiffFormulaeMod_gr_186.gif],  and an even number of  [Graphics:Images/NumericalDiffFormulaeMod_gr_187.gif] equations.  Since we want to include the remainder terms, we need to use series expansions of order  [Graphics:Images/NumericalDiffFormulaeMod_gr_188.gif].  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 [Graphics:Images/NumericalDiffFormulaeMod_gr_189.gif] instead of  [Graphics:Images/NumericalDiffFormulaeMod_gr_190.gif].  

[Graphics:Images/NumericalDiffFormulaeMod_gr_191.gif]

The Five Point Central Difference Formulae

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

[Graphics:Images/NumericalDiffFormulaeMod_gr_192.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_193.gif]

 

Higher order formulas are easy to obtain with Mathematica.
Exploration

 

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.  

[Graphics:Images/NumericalDiffFormulaeMod_gr_198.gif]
[Graphics:Images/NumericalDiffFormulaeMod_gr_199.gif]

Explorations with the Higher Order Formulas

[Graphics:Images/NumericalDiffFormulaeMod_gr_200.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_201.gif]

 

[Graphics:Images/NumericalDiffFormulaeMod_gr_202.gif]

[Graphics:Images/NumericalDiffFormulaeMod_gr_203.gif]

 

Higher order formulas are easy to obtain with Mathematica.
Exploration

 

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.  

 

Download this Mathematica Notebook Numerical Differentiation Formulae

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(c) John H. Mathews 2004