This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: gsl_root_test_residual(double f, double epsabs)


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).




Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]