|
From: | Fotios |
Subject: | Re: "Exact" jacobian calculation |
Date: | Sun, 02 Oct 2011 12:24:52 +0300 |
User-agent: | Mozilla/5.0 (Windows NT 6.1; WOW64; rv:7.0.1) Gecko/20110929 Thunderbird/7.0.1 |
Στις 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))/hHope 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
[Prev in Thread] | Current Thread | [Next in Thread] |