[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [Bug-gsl] diff.c (numerical differentiation)

**From**: |
spasmous |

**Subject**: |
[Bug-gsl] diff.c (numerical differentiation) |

**Date**: |
Sun, 6 Apr 2008 13:29:41 -0700 |

In several places in the diff.c code there are lines like the following:
h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
...
*result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
There is an assumption here that the difference between (x+h) and x is
exactly identical to h, which is true only in infinite precision. In finite
precision there can be a small error introduced because x+h may not be
exactly representable in finite precision and so a nearby number is used
(rounded or truncated). It is very easy to remove this error in the
following way:
h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
...
h = (x+h)-x; // ensure difference between (x+h) and x is exactly h
*result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
The issue is explained in detail here -
http://www.nrbook.com/a/bookcpdf/c5-7.pdf

**[Bug-gsl] diff.c (numerical differentiation)**,
*spasmous* **<=**