Re: Numerical differentiation, recursion of functions not working

[Martin Helm]

Am Montag, 29. März 2010 14:42:37 schrieb James Sherman Jr.:
> If you want to do it more "mathematically" you should code up and actual
> derivative operator/function which takes a function as its input and a
> function as its output.  (Or since you want to do it numerically, a set of
> points instead of a function.)  Like Thomas said, your dx function is "the
> derivative of f evaluated at x" so when you compose this function with
> itself you get "the derivative of f evaluated at (the derivative of f
> evaluated at x)"
> 
> So, coding wise it would look like (assuming you call the derivative
> function "deriv")
> 
> fprime = deriv(f);
> fprime2 = deriv(fprime)
> 
> Then, fprime2(x) is what you're looking for.
> 
> On Mon, Mar 29, 2010 at 6:32 AM, Stefan Neumann <address@hidden> wrote:
> > 2010/3/29 Thomas Shores <address@hidden>
> >
> >  The problem is in your coding.  You have f(x) hardwired into dx, so
> > you're
> >
> >> not really recursing on a finite difference formula when you call it via
> >> dx(dx(x)).  What you're actually calculating is the first difference of
> >> f(x) evaluated at the first difference of f(x).
> >
> > Yes, I noticed.
> >
> > Coding-wise dx(dx(x)) does not work.
> >
> > What bugs me is that mathematically the second derivative does equal the
> > first derivative of the first derivative.
> > So it should be possible to code that using a language modeled on
> > mathematical thinking.
> >
> > Obviously you can always work around the problem by using the established
> > formulas with differential quotients, but you have to use different
> > formulas for 1st, 2nd, 3rd etc derivative.
> >
> > Not really elegant.
> > So I was wondering how an elegant solution would look.
> >
> > Regards,
> > Stefan
> >

You can (just for fun) do some lisp like higher order function trickery, for 
example

function df = deriv(f)
  df = @(x) (f(x+1e-5)-f(x-1e-5))/2e-5;
endfunction

df = deriv(@sin)

=> df =
@(x) (f (x + 1e-5) - f (x - 1e-5)) / 2e-5


d2f = deriv(df)

=> d2f =
@(x) (f (x + 1e-5) - f (x - 1e-5)) / 2e-5

This seems to look the same but is not   :-)

df(0)
=> 1.00000

d2f(0)
=> 0

as it should be!

Think about it (the trick is I am using closures here by using anonymous 
functions - or more lispy: lambda expressions)

- mh