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 <awsteiner@.........> 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;
>
