This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
bug in root finding?
- To: gsl-discuss at sources dot redhat dot com
- Subject: bug in root finding?
- From: Jon Danielsson <jond at rhi dot hi dot is>
- Date: Tue, 13 Nov 2001 16:44:17 +0000 (GMT)
I am using gsl 1.0 on RH7.2, and encountered the problem that when
the rootfinder (gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
finds the solution immeately, i.e. its brackeded between -100,100
on third function evaluaiton it zero, which is the exact solution, but
then fails to recognize the fact, loops until max_iter is reached
here is the codefragment
gsl_root_fsolver *s;
gsl_function F;
struct FRho_params params = {b,lambda,rbar,gamma};
F.function = &FRho;
F.params = ¶ms;
params.b=b;// = {d1,lambda,rbar,gamma};
s = gsl_root_fsolver_alloc (gsl_root_fsolver_brent);
gsl_root_fsolver_set (s,&F, -100,100);
do {
iterations++;
status = gsl_root_fsolver_iterate (s);
r = 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.00001);
}
while (status == GSL_CONTINUE && iterations < max_iterations);
and the function
nnuble FRho (double x, void *params) {
struct FRho_params *p = (struct FRho_params *) params;
double b= p->b;
double lambda = p->lambda;
double rbar= p->rbar;
double gamma= p->gamma;
double res=0.0;
res= (x+lambda*(1-b)) - lambda * FN(sqrt(gamma)*(x-rbar));
return res;
}
double FN(double x){
return 0.5+ gsl_sf_erf (x/sqrt(2.0))/(2.0);
}