[Bug-gsl] rk8pd

[Andrew W. Steiner]

I'm a bit surprised by the lack of accuracy of the 8th order stepper in a
simple case, and I'm curious if anyone else has run into this for
other problems.
Admittedly, "higher order does not always mean higher accuracy", but if this
behavior is consistent over several functions it points to a possible bug.
In any case, in the code below I use the example
from the GSL manual and applies it to a Bessel function and compare the
PD and CK steppers. I consistently find the latter outperforming the
former, usually
by a factor of 2.

Any ideas?

Thanks,
Andrew

int gsl_derivs(double x, const double y[], double dydx[], void *vp) {
  dydx[0]=y[1];
  if (x==0.0) dydx[1]=0.0;
  else dydx[1]=(-x*y[1]+(-x*x+1)*y[0])/x/x;
  return 0;
}

  {
    const gsl_odeiv_step_type * T  = gsl_odeiv_step_rkck;
  // const gsl_odeiv_step_type * T  = gsl_odeiv_step_rk8pd;

    gsl_odeiv_step * s     = gsl_odeiv_step_alloc (T, 2);

    gsl_odeiv_system sys = {gsl_derivs, 0, 2, 0};

    double tt = 0.0, t1 = 1.0;
    double h = 1e-1;
    double ya[2]={0.0,0.5};
    double y_err[2];
    double dydt_in[2], dydt_out[2];

    /* initialise dydt_in from system parameters */
    GSL_ODEIV_FN_EVAL(&sys, tt, ya, dydt_in);

    while (tt < t1) {
      int status = gsl_odeiv_step_apply
        (s, tt, h, ya, y_err, dydt_in, dydt_out, &sys);
      if (status != GSL_SUCCESS) break;
      dydt_in[0] = dydt_out[0];
      dydt_in[1] = dydt_out[1];
      tt += h;
      //printf ("%.5e %.5e %.5e\n", tt, ya[0], ya[1]);
      cout << tt << " " << ya[0] << " " << ya[1] << " ";
      cout << fabs((ya[0]-gsl_sf_bessel_J1(tt))/gsl_sf_bessel_J1(tt)) << endl;
    }

    gsl_odeiv_step_free (s);
  }
  cout << endl;