Re: [Help-gsl] Jacobian evaluation in ODE suite

[Dale Lukas Peterson]

On Thu, Sep 23, 2010 at 05:05:59PM +0200, Claude Lacoursiere wrote:
> 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

Claude,
  I tried your example and found the same thing.  I am also surprised
  by this behavior.  The manual is vague about which of the
  solvers make use of the Jacobian.  It seems to me that only the
  implicit methods should use it, but maybe this isn't correct?  Even
  if this is the case, the methods you used should make use of the
  Jacobian when they employ Newton's method to solve the system of
  equations, right?

  I had a look at the source files and it seems that G. Jungman is the
  author of most (if not all) of the ode-initval code.  Maybe we could
  contact him directly if he doesn't respond on here?

Cheers,
Luke

> #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;
> }
> 
>