[Bug-gsl] Re: rk8pd

[Andrew W. Steiner]

Ok. I checked a couple other functions and they're fine.
 It seems I just picked a case where the higher-order
stepper does worse. There's probably some deep reason why the Bessel
function is
difficult at higher order. I have a couple guesses, but if someone
wants to enlighten
me that'd be great too.

On Mon, Apr 7, 2008 at 9:09 AM, Andrew W. Steiner <address@hidden> wrote:
> 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;
>