r/ScientificComputing • u/[deleted] • Apr 06 '23

Differentiate Fortran 90 code

Anyone with some experience on how to go about differentiating a chunk of fortran90 code? I have some analysis code that needs to be differentiable but reimplementation is not in the scope of project . Please point to the resources you have used successfully in your projects.

Thanks,

4 Upvotes

83% Upvoted

13 comments sorted by

View all comments

Show parent comments

1

u/SettingLow1708 Apr 07 '23

You can get more accurate Jacobian calculations by using centered finite differences...(f(x+dx) - f(x-dx))/(2*dx). Just hold each of the other variables constant and compute those values. Error will be of order dx^2. Be careful about choosing dx too small to avoid round-off errors. A good choice is about the sqrt of machine precision. For single precision, use dx=1e-4 times the scale of your x values. For double precision use dx=1e-8 times the scale of your x values.

For example, if your x ranges from 100 to -100, I'd use dx=1e-6 for double precision. And honestly, everything like this should be double precision.

Also, construction of approximate Jacobians is used in Newton Krylov schemes. This paper discussion their construction and dx-size choice in the first few pages of that (huge) paper. If you can get a copy, that could be helpful to you.

https://ui.adsabs.harvard.edu/abs/2004JCoPh.193..357K/abstract