This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
gsl_diff_central
- From: Andrew Steiner <stein at physics dot umn dot edu>
- To: gsl-discuss at sources dot redhat dot com
- Date: Thu, 19 Jun 2003 09:44:47 -0500 (CDT)
- Subject: gsl_diff_central
Hello!
I had some ideas for changes to gsl_diff_central in response to my
earlier complaint that it fails for large values of "x" (BUG #13). I
wanted to post them so that someone could check and make sure that I
haven't missed something. If people like it, then I'll make a proper
patch or whatever.
int gsl_diff_central(const gsl_function *f, double x,
double *result, double *abserr) {
int i, k;
// Change the stepsize a little
double h=GSL_SQRT_DBL_EPSILON*x;
double a[4],d[4],a3;
for (i=0;i<4;i++) {
a[i]=x+(i-2.0)*h;
d[i]=GSL_FN_EVAL(f,a[i]);
}
// Changed k<5 to k<4
for (k=1;k<4;k++) {
for(i=0;i<4-k;i++) {
d[i]=(d[i+1]-d[i])/(a[i+k]-a[i]);
}
}
// Here I changed d[0]+d[1]+d[2]+d[3] to d[0].
// d[0] contains the third derivative and the other entries
// may have different units
a3=fabs(d[0]);
// Calculate the new stepsize
h=pow(GSL_SQRT_DBL_EPSILON*x/(2.0*a3),1.0/3.0);
if (h>100.0*GSL_SQRT_DBL_EPSILON*x) {
h=100.0*GSL_SQRT_DBL_EPSILON*x;
}
*result=(GSL_FN_EVAL(f,x+h)-GSL_FN_EVAL(f,x-h))/(2.0*h);
*abserr=fabs(100.0*a3*h*h);
// We can't restrict a3 directly since we don't know the units
// for f'''(x), so we ensure that the error is not too much
// smaller than the result instead.
if ((*abserr)<(*result)*100.0*GSL_SQRT_DBL_EPSILON) {
(*abserr)=(*result)*100.0*GSL_SQRT_DBL_EPSILON;
}
return GSL_SUCCESS;
}
Thanks,
Andrew
----------------------------------------------------------------------
Andrew W. Steiner Post-doctoral Research Associate
Nuclear Theory Group University of Minnesota
Phone: 612-624-7872 Fax: 612-624-4875
Email: stein@physics.umn.edu URL: http://umn.edu/~stein178
----------------------------------------------------------------------