Friday, May 20, 2011

Numerical Differentiation using Mathematica

I am teaching myself Mathematica over summer, and coded up a simple module that computes numerical differentiation formulas automatically. Essentially, the module enables one to compute arbitrary differentiation rules like those explained here.

This module spits out the different formulae for a m-point approximation to the n-th derivative, and the leading error term.

Here is the (updated) module:

AppDeriv[m_, n_] := Module[{points, IP, nthDeriv, e1},
  points = Table[{x0 + i h, Subscript[f, i]}, {i, 0, m - 1}];
  IP = InterpolatingPolynomial[points, x];
  nthDeriv = D[IP, {x, n}];
  e1 = (1/Factorial[m])*
    D[Apply[Times, (x - Table[{x0 + i h}, {i, 0, m - 1}])], {x, n}];
  Formula =
    Table[{Subsuperscript[f, i, n], nthDeriv /. x -> (x0 + i h),
       Superscript[f, m] e1 /. x -> (x0 + i h)}, {i, 0, m - 1}] //
     Simplify ]]

This seems to work alright for me (click to enlarge).


Anonymous said...

AppDeriv[2, 1]

InterpolatingPolynomial::list: List expected at position {At Line = 3, the input was:,AppDeriv[2,1],1} in {At Line = 3, the input was:,AppDeriv[2,1],InterpolatingPolynomial[points,x]}. >>

You might want to learn programming a little more before posting such code.

Sachin Shanbhag said...

You're right. I think I fixed it now. As bonus, it also gives the leading error term.