[Help-gsl] solving the Van der Pol problem with rk1imp with fixed step size using low level functions

[Juan Pablo Amorocho D.]

Hi all,

I'm trying to solve the Van der Pol problem that is in the gsl documentation
using the step type  (aka solver) rk1imp - Euler backward - with fixed step
size and low level functions (evolve, step, control). So far all I get is
that the status from the evolve_apply_fixed_step is 3. Looking that the gsl
code I found that this means that there is pointer problem. Now, if I change
the solver to, say, rk2 the code runs. The strange thing is that if I use
the driver with fixed step and rk1imp, the code runs.  Any ideas of what may
be wrong? Any help would be appreciated.

-- Juan

PS. Here is my code using the low level functions:

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

 int func (double t, const double y[], double f[],
           void *params)
     {
       double mu = *(double *)params;
       f[0] = y[1];
       f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
       return GSL_SUCCESS;
     }

     int jacobian (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_matrix * m = &dfdy_mat.matrix;
       gsl_matrix_set (m, 0, 0, 0.0);
       gsl_matrix_set (m, 0, 1, 1.0);
       gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
       gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
       dfdt[0] = 0.0;
       dfdt[1] = 0.0;
       return GSL_SUCCESS;
     }


  int main (void)
  {
    const double epsAbs = 1e-4, epsRel = 1e-4;

    double mu = 10.;
    double y[2] = { 1.0, 0.0 };

    double t = 0.0, tf = 10.0, h = 1e-4, hstart = h;
    int status = 0;

    gsl_odeiv2_system system = { func, jacobian, 2, &mu };

    const gsl_odeiv2_step_type * solvingAlgorithm = gsl_odeiv2_step_rk1imp;

    gsl_odeiv2_step * step = gsl_odeiv2_step_alloc(solvingAlgorithm, 2);

    gsl_odeiv2_control * control = gsl_odeiv2_control_y_new(epsAbs, epsRel);

    gsl_odeiv2_evolve * evolve = gsl_odeiv2_evolve_alloc(2);


    while (t < tf)
      {
        status = gsl_odeiv2_evolve_apply_fixed_step(evolve, control, step,
&system, &t,hstart, y);

        if (status != GSL_SUCCESS)
          {
            printf ("error: evolve returned %d\n", status);
            break;
          }

        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
      }

    gsl_odeiv2_evolve_free(evolve);
    gsl_odeiv2_control_free(control);
    gsl_odeiv2_step_free(step);
    return 0;
  }


PSS And here my code using the driver:

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

 int func (double t, const double y[], double f[],
           void *params)
     {
       double mu = *(double *)params;
       f[0] = y[1];
       f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
       return GSL_SUCCESS;
     }

     int jac (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_matrix * m = &dfdy_mat.matrix;
       gsl_matrix_set (m, 0, 0, 0.0);
       gsl_matrix_set (m, 0, 1, 1.0);
       gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
       gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
       dfdt[0] = 0.0;
       dfdt[1] = 0.0;
       return GSL_SUCCESS;
     }


  int main (void)
  {
    const double hstart = 1e-5, epsAbs = 1e-4, epsRel = 1e-4;
    double mu = 10.;
    const unsigned long numberOfSteps = 10000;
    const int problemDimension = 2;
    gsl_odeiv2_system sys = { func, jac, problemDimension, &mu };

    gsl_odeiv2_driver *d =
      gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk1imp,
                                     hstart, epsAbs, epsRel);

    double t = 0.0;
    double y[2] = { 1.0, 0.0 };
    int i, s;

    for (i = 0; i < numberOfSteps; i++)
      {
        s = gsl_odeiv2_driver_apply_fixed_step (d, &t, hstart,
numberOfSteps, y);

        if (s != GSL_SUCCESS)
          {
            printf ("error: driver returned %d\n", s);
            break;
          }

        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
      }

    gsl_odeiv2_driver_free (d);
    return s;
  }