Re: "Exact" jacobian calculation

[Fotios]

Στις 2011-10-02 10:46, ο/η Søren Hauberg έγραψε:
> søn, 02 10 2011 kl. 06:05 +0300, skrev Fotios:
> > I attach a simple function with the hope that someone will find it
> > useful. It calculates the "exact" jacobian of an analytic function
> > f:R^n->R^n avoiding subtractive cancellation by following the imaginary
> > step approach. Any corrections are welcomed. Enjoy.
> This looks interesting (you should consider finding a place for it in
> Octave-Forge). Two comments:
>
> * I think you should pre-allocate 'Df' before the loop.
> * Instead of writing 'Df = Df  / h' you can write 'Df /= h', which might
> be faster.
>
> I have never heard of the Complex Step Method, but it seems similar to
> finite differences. Does it make sense have different magnitudes for
> different coordinates (i.e. to let 'h' be a vector rather than a single
> number) ? (I'm asking out of ignorance)
>
> Søren

Indeed the dfdx should be Df i just changed my notation and i forgot 
that one. The advantage of the method compared to fd is that since there 
is no subtractive cancellation you can essentially choose your step to 
be as small as you want it, therefore the default value is h=1e-20. I do 
not think it makes sense to choose h to be vector for a generic function 
like that and for so small h value, you always (?) get the "exact" value 
up to machine epsilon. Regarding the method for a function f:R->R you have

f(x+ih) = f(x) + ihf'(x) + (ih)^2 f''(x) / 2 + ... => [take imaginary part]
im (f(x+ih)) = im (f(x) + ihf'(x) + (ih)^2 f''(x) / 2 + ...) =>
im(f(x+ih)) = hf'(x) [+ O(h^3)] =>
f'(x) \approx im(f(x+ih))/h

Hope it helps a bit. This function and some others will be part of a pkg 
for chaotic dynamics that i prepare. Thanks for commenting Søren.

/Fotis