This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl_root_test_residual(double f, double epsabs)
- From: Karsten Howes <karsten at videotron dot ca>
- To: Brian Gough <bjg at network-theory dot co dot uk>
- Cc: gsl-discuss at sources dot redhat dot com
- Date: Thu, 09 Sep 2004 08:37:17 -0400
- Subject: Re: gsl_root_test_residual(double f, double epsabs)
- References: <40C3AF75.7080904@yahoo.es> <40C4229F.3070907@yahoo.es><opsdq60tes0c67l8@localhost> <16703.16611.701140.367487@network-theory.co.uk>
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).