ODE solver

[Ivo Alxneit <ivo dot alxneit at psi dot ch>]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

dear everybody

i tried to use the various ODE solvers provided by gsl (version 1.4) and i got 
a bit confused. according to the manual "some" solvers make use of the 
jacobian. it seems that "some" until version 1.4 is "bsimp only" (maybe it 
would help to explicitly state in the manual if some algorithm does not make 
use of the jacobian).

to convince myself it tested all algorithms provided with the demo provided in 
the manual (van der pol, gsl_odeiv_evolve). i tried two cases for the 
jacobian (jac_with corresponds to the manual).

int
jac_with (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_vector_view dfdt_vec = gsl_vector_view_array (dfdt, 2);
  gsl_matrix_set (&dfdy_mat.matrix, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat.matrix, 0, 1, 1.0);
  gsl_matrix_set (&dfdy_mat.matrix, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (&dfdy_mat.matrix, 1, 1, -mu*(y[0]*y[0] - 1.0));
  gsl_vector_set (&dfdt_vec.vector, 0, 0.0);
  gsl_vector_set (&dfdt_vec.vector, 1, 0.0);
  return GSL_SUCCESS;
}
int
jac_without (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  dfdy = NULL;
  dfdt = NULL;
  return GSL_SUCCESS;
}

is jac_without correctly used?

then again i got the following results (the functions look correct in ALL 
cases) for the number of steps needed:

Algorithm: rk2
number of steps (with jacobian): 11714
number of steps (without jacobian): 11714

Algorithm: rk4
number of steps (with jacobian): 84197
number of steps (without jacobian): 84197

Algorithm: rkf45
number of steps (with jacobian): 1495
number of steps (without jacobian): 1495

Algorithm: rkck
number of steps (with jacobian): 1215
number of steps (without jacobian): 1215

Algorithm: rk8pd
number of steps (with jacobian): 722
number of steps (without jacobian): 722

Algorithm: rk2imp
number of steps (with jacobian): 82645
number of steps (without jacobian): 82645

Algorithm: rk4imp
number of steps (with jacobian): 84270
number of steps (without jacobian): 84270

Algorithm: bsimp
number of steps (with jacobian): 144
number of steps (without jacobian): 262

Algorithm: gear1
number of steps (with jacobian): 82570
number of steps (without jacobian): 82570

Algorithm: gear2
number of steps (with jacobian): 72503
number of steps (without jacobian): 72503

so what happens with bsimp if the jacobian is not provided (jac_without). 
according to the manual i would have expected that i get an error message or 
that the thing bombs out due to the null pointer.
- -- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
Paul Scherrer Institute             fax: +41 56 310 2688
CH-5232 Villigen                   http://solar.web.psi.ch
Switzerland                   gnupg key: 0x515E30C7
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.3 (GNU/Linux)

iD8DBQFBf3INAd7CE1FeMMcRAl2gAKCWWOdT5rTkx2xxsZhClXmTAhzLLQCgpu4m
8qpPqj+qzz7J5IQf8KMalDQ=
=17ZC
-----END PGP SIGNATURE-----