[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] bug report regarding gsl_diff_*

From: Zegenhagen Stefan
Subject: [Bug-gsl] bug report regarding gsl_diff_*
Date: 23 Feb 2004 09:02:41 +0100

Dear all,

I think I've found a bug that makes the numerical differentiation
routines quite unusable in some cases.

As visible from the source file diff/diff.c, in all the differentiation
routines, you calculate a higher order derivative first to choose a
stepsize that minimizes the absolute error when calculating f'. As a
stepsize to calculate the higher order derivative, you always choose the
absolute value of


This method fails if the value x at which to calculate the derivative is
so large that (x + i*GSL_SQRT_DBL_EPSILON) == x holds due to the
restricted numerical precision. Then, some of the intermediate values in
the calculation will become infinite (because of division by zero) and
the final result for h is NaN in the calculation of f'. This is

I also think that the same holds if x is so small that the equation (x +
i*GSL_SQRT_DBL_EPSILON) == GSL_SQRT_DBL_EPSILON holds. Both conditions
may be fulfilled even for double precision values.

The solution could be to either a stepsize relative to x in the
calculation of the higher order derivative(something like 

        double h = x * GSL_SQRT_DBL_EPSILON;

in the function header). This solution should avoid the restricted
precision problems.

Sincerely yours,

Dipl. Phys. Stefan Zegenhagen
Max-Planck-Institut fuer Plasmaphysik
- Teilinstitut Greifswald -
Wendelsteinstr. 1
D-17491 Greifswald, Germany
Tel. +49 (0)3834 88 1277   Fax:  +49 (0)3834 88 2509
e-mail: address@hidden

reply via email to

[Prev in Thread] Current Thread [Next in Thread]