Hello.
I'm using the ODE module in a class I teach and I was surprised that the
Jacobian function is unused for all integration methods with the exception of
gsl_odeiv_step_bsimp
I tried all methods on the BZ reaction which is very stiff, and none of
gsl_odeiv_step_rk2imp
gsl_odeiv_step_rk4imp
gsl_odeiv_step_gear2
used the provided Jacobian function a single time.
Using either the Matlab toolbox or the Octave odepkg, all implicit methods need
to evaluate the Jacobian to produce a decent solution for this
equation. True, neither
of these have the Gauss' colocation method available in GSL.
Does anyone know how to force the evaluation of Jacobians? Is this a bug?
Attached is the code if anyone is curious.
Best,
/Claude
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
typedef struct bz_data {
int count;
int jac_count;
double alpha;
double beta;
double gamma;
} bz_data;
int bz (double t, const double y[], double f[], void *params){
bz_data * b = (bz_data *) params;
f[0] = b->alpha*(y[1] + y[0]*(1-b->beta*y[0] - y[1]));
f[1] = (1/b->alpha)*(y[2] - (1+y[0])*y[1]);
f[2] = b->gamma*(y[0] - y[2]);
++b->count;
return GSL_SUCCESS;
}
int jac_bz (double t, const double y[], double *dfdy, double dfdt[], void *params) {
bz_data *p = (bz_data *)params;
double alpha = p->alpha;
double beta = p->beta;
double gamma = p->gamma;
gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 3, 3);
gsl_matrix * m = &dfdy_mat.matrix;
p->jac_count ++ ;
gsl_matrix_set (m, 0, 0, alpha*(1-beta*y[0]-y[1] - beta*y[0]));
gsl_matrix_set (m, 0, 1, alpha*(1 - y[0]));
gsl_matrix_set (m, 0, 2, 0);
gsl_matrix_set (m, 1, 0, (1/alpha)*(-y[1]));
gsl_matrix_set (m, 1, 1, -(1/alpha)*(1+y[0]));
gsl_matrix_set (m, 1, 2, (1/alpha));
gsl_matrix_set (m, 2, 0, gamma*y[0]);
gsl_matrix_set (m, 2, 1, 0);
gsl_matrix_set (m, 2, 2, -gamma*y[2]);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
dfdt[2] = 0.0;
return GSL_SUCCESS;
}
int main () {
//const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2imp;
const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4imp;
//const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
//const gsl_odeiv_step_type * T = gsl_odeiv_step_bsimp;
//const gsl_odeiv_step_type * T = gsl_odeiv_step_gear2;
gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 3);
gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (3);
bz_data dat = {0, 0, 77.27, 8.375E-6, 1.161};
gsl_odeiv_system sys = {bz, jac_bz, 3, &dat};
double t = 0.0, t1 = 360;
double h = 1e-6;
double y[] = { 1, 2, 3};
double t2 = t;
double interval = 0.05; // suppress excessive output
while (t < t1)
{
int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
if (status != GSL_SUCCESS)
break;
if (t > t2 + interval ) {
printf ("%.5e %.5e %.5e %5e %d\n", t, y[0], y[1], y[2], dat.count);
t2 = t;
}
}
gsl_odeiv_evolve_free (e);
gsl_odeiv_control_free (c);
gsl_odeiv_step_free (s);
printf("Number of Jacobian evaluations = %d\n"
"Number of Function evaluations = %d\n", dat.jac_count,
dat.count);
return 0;
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
