gsl_root_test_residual(double f, double epsabs)

[Karsten Howes]

OK, thanks Brian.  I am using Brent's method which appears always to have 
valued the function at the root.  I solved the problem by storing a map 
root->value inside my function.  I have included the code below.  It's a 
C++/stl/boost wrapper for the brent solver, which will be of no interest 
to those who are hard core C programmers. And, it has a few pointers of 
overhead, etc, which doesnt matter to me in my applications.  It does have 
the advantage of being able to find roots of member functions though 
(using e.g. boost::bind).

Best regards,
Karsten

   //params for gsl root finder.  Includes the function
   struct DoubleGslBoostFuncWrapper {
     boost::function<double (double)> f;
     std::map<double, double> f_map;
   };

   static double doubleGslFunc(double x, void *boostFunc){
     double f = ((DoubleGslBoostFuncWrapper *) boostFunc)->f(x);
     ((DoubleGslBoostFuncWrapper *) boostFunc)->f_map[x] = f;
     return f;
   }

   //Finds root of f with error of tolerance
   int brent(boost::function<double (double)> &f,
	    double x_lo, double x_hi, double tolerance, double &root){
     int status;
     int iter = 0, max_iter = 20;
     const gsl_root_fsolver_type *T;
     gsl_root_fsolver *s;
     gsl_function F;

     DoubleGslBoostFuncWrapper gslBoostFuncWrapper;
     gslBoostFuncWrapper.f = f;
     F.function = &doubleGslFunc;
     F.params = &gslBoostFuncWrapper;

     T = gsl_root_fsolver_brent;
     s = gsl_root_fsolver_alloc (T);
     gsl_root_fsolver_set (s, &F, x_lo, x_hi);

     printf ("using %s method\n",
	    gsl_root_fsolver_name (s));

     printf ("%5s [%9s, %9s] %9s %10s %9s\n",
	    "iter", "lower", "upper", "root",
	    "err", "err(est)");

     do
       {
	iter++;
	status = gsl_root_fsolver_iterate (s);
	root = gsl_root_fsolver_root (s);
	x_lo = gsl_root_fsolver_x_lower (s);
	x_hi = gsl_root_fsolver_x_upper (s);
	//status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001);
	double residual = gslBoostFuncWrapper.f_map[root];
	status = gsl_root_test_residual (residual, tolerance);

	if (status == GSL_SUCCESS)
	  printf ("Converged:\n");

	printf ("%5d [%.7f, %.7f] %.7f %.7f, f=%.7f\n",
		iter, x_lo, x_hi, root, x_hi - x_lo, residual);
       }
     while (status == GSL_CONTINUE && iter < max_iter);
     return status;
   }



On Wed, 08 Sep 2004 18:26:59 +0100, Brian Gough <bjg@network-theory.co.uk> 
wrote:

> Karsten Howes writes:
>  > I am working with the 1D root finder framework and I am trying to 
> figure
>  > out how to get the value of my function at the current best root 
> estimate
>  > (without re-evaluating my function of course).
>
> You do need to evaluate the function f(root).  Depending on the
> algorithm it may not be a reevaluation (see e.g. bisection).
>