![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhkT04jxOOu6D_83ndvEh20R-GXBuhqUxobkxF46WCxuWUxqmrGzIYbLPM4_ThaDZ2oZ2giLB7yFa2GMwNVo6OU5x_dwpCd2o2w3VU_xDhIw8q8x75bMq5xAkm1g6hcFpehVMtZU0gU9ek/s400/2d-gaussian.png)
A finite difference approximation would be (in Octave):
dfdx = ( f(x+h) - f(x) )/ h;
Complex step is even cooler, because you don't have any subtraction (there's no difference), so you can choose a very small step size without loosing accuracy due to subtractive cancelation:
complex_step = complex( 0, 1e-320 );
The derivative is approximated by just the imaginary part:
x = x + complex_step;
dfdx = imag( f(x) )/imag( complex_step );
So the derivative with respect to x looks like this:
![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBg2n9ULJ3e2KOB5TG6uheJYFJgvXExTrHJtfX9lmGeBLS3LQlV_n41n_mnjPLlSePnPAHrX6jr2Dp9IEbXUJH5DvO2zWn7Ljcg4HFSbQMNYrRu_IqxnVFerp4k874vccKvcVQr65dtdo/s400/2-d_gauss_dy.png)
And the estimated derivative with respect to y:
![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEha5lD3g7T-1FHaLK5x-Er1Q-OrZmVczgGbHFmQAbX2R2dC9jwlNr2O98AdPggWyZZ0CK6awJ4CIo0FrXnrJwAcZRT3SCabiscrH-KeX-iFN7MM3mW2vi-n-eji6QWjsmduRQDUnAnzd-8/s400/2-d_gauss_dx.png)
This approach is really useful for design sensitivity analysis, and since modern Fortran supports complex types we can even use this method for serious number crunching!
No comments:
Post a Comment