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 _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl