gsl_diff_central

[Andrew Steiner <stein at physics dot umn dot edu>]

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