This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
ODE solver
- From: Ivo Alxneit <ivo dot alxneit at psi dot ch>
- To: gsl-discuss at sources dot redhat dot com
- Date: Wed, 27 Oct 2004 12:01:40 +0200
- Subject: ODE solver
- Organization: Paul Scherrer Institute
-----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-----