[Help-gsl] Bug (or not?) in multiroots/fdjacobian.c

[eknecronzontas]

Since I can't quite decide if this is a bug, I'll post
it here first...

The stepsize for finite-differencing in
gsl-1.8/multiroots/fdjacobian.c is specified by the
lines

> double xj = gsl_vector_get (x, j);
> double dx = epsrel * fabs (xj);

This is, of course mathematically correct.
Nevertheless, the behavior is less than ideal if one
of the elements of the vector 'x' is sufficiently
small. This occurs often if one component of the root
happens to be zero. If this occurs, then 'dx' can be
so small that one of the columns of the Jacobian is
identically zero. This immediately leads to NAN's
which are not easy to trace.

The following code:
> if (fddx == 0) {
>   fddx = epsrel;
> }
does not fix the issue, as fddx can be finite and very
small, yet the change in the equations, (g1-g0)/fddx,
is zero.

This will cause a Jacobian with zero determinant even
for a trival function, so long as the initial guess is
sufficiently small, e.g. 1.0e-20.

Thus I'm convinced that fdjacobian.c ought to be
modified. One naive guess is 

> double dx = epsrel * fabs (xj) + epsabs;

Alternatively, one could just add the line

> if (dx < epsmin) dx = epsmin;

somewhere later in fdjacobian(). However, both of
these solutions would introduce a scale to the
problem, which might not be good. Argh.

Alternatively, most of the difficulty might be settled
by simply calling the error handler in set()-like
functions if one of the columns in the Jacobian is all
zeros. This would at least warn the user that
iterate() will fail.

Finally, a last option (most difficult) would be to
allow the user to specify the stepsize to be taken for
each component of the initial guess vector. 

Let me know what you think,
Andrew Steiner