GSL ODE Solver Example

ode_test.c

/*
  Example adapted from the GNU Scientific Library Reference Manual
  Edition 1.1, for GSL Version 1.1
  9 January 2002
  URL: gsl/ref/gsl-ref_25.html#SEC381
   
  Revisions by:  Dick Furnstahl  [email protected]
                 Pedro Rittner   [email protected]
 
  Revision history:
  11/03/02  changes to original version from GSL manual
  11/14/02  added more comments
  06/09/14  changed build suggestions to comply with modern gcc linkage for GSL/BLAS
*/

/* 
   Compile and link with:
   gcc -c ode_test.c
   gcc -o ode_test ode_test.o -lgsl -lgslcblas -lm

   Alternatively:
   gcc -std=c99 -Wall -O0 -g3 ode_test.c -o ode_test.bin -lgsl -lblas
*/

/* 
   The following program solves the second-order nonlinear 
   Van der Pol oscillator equation (see Landau/Paez 14.12, part 1),

   x"(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0

   This can be converted into a first order system suitable for 
   use with the library by introducing a separate variable for 
   the velocity, v = x'(t).  We assign x --> y[0] and v --> y[1].
   So the equations are:
   x' = v                  ==>  dy[0]/dt = f[0] = y[1]
   v' = -x + \mu v (1-x^2) ==>  dy[1]/dt = f[1] = -y[0] + mu*y[1]*(1-y[0]*y[0])
*/

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

/* function prototypes */
int rhs (double t, const double y[], double f[], void *params_ptr);
int jacobian (double t, const double y[], double *dfdy,
              double dfdt[], void *params_ptr);

/*************************** main program ****************************/

int
main (void)
{
    int dimension = 2;		/* number of differential equations */
  
    double eps_abs = 1.e-8;	/* absolute error requested */
    double eps_rel = 1.e-10;	/* relative error requested */

    /* define the type of routine for making steps: */
    const gsl_odeiv_step_type *type_ptr = gsl_odeiv_step_rkf45;
    /* some other possibilities (see GSL manual):          
       = gsl_odeiv_step_rk4;
       = gsl_odeiv_step_rkck;
       = gsl_odeiv_step_rk8pd;
       = gsl_odeiv_step_rk4imp;
       = gsl_odeiv_step_bsimp;  
       = gsl_odeiv_step_gear1;
       = gsl_odeiv_step_gear2;
    */

    /* 
       allocate/initialize the stepper, the control function, and the
       evolution function.
    */
    gsl_odeiv_step *step_ptr = gsl_odeiv_step_alloc (type_ptr, dimension);
    gsl_odeiv_control *control_ptr = gsl_odeiv_control_y_new (eps_abs, eps_rel);
    gsl_odeiv_evolve *evolve_ptr = gsl_odeiv_evolve_alloc (dimension);

    gsl_odeiv_system my_system;	/* structure with the rhs function, etc. */

    double mu = 10;		/* parameter for the diffeq */
    double y[2];			/* current solution vector */

    double t, t_next;		/* current and next independent variable */
    double tmin, tmax, delta_t;	/* range of t and step size for output */

    double h = 1e-6;		/* starting step size for ode solver */

    /* load values into the my_system structure */
    my_system.function = rhs;	/* the right-hand-side functions dy[i]/dt */
    my_system.jacobian = jacobian;	/* the Jacobian df[i]/dy[j] */
    my_system.dimension = dimension;	/* number of diffeq's */
    my_system.params = &mu;	/* parameters to pass to rhs and jacobian */

    tmin = 0.;			/* starting t value */
    tmax = 100.;			/* final t value */
    delta_t = 1.;

    y[0] = 1.;			/* initial x value */
    y[1] = 0.;			/* initial v value */

    t = tmin;             /* initialize t */
    printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);	/* initial values */

    /* step to tmax from tmin */
    for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)
    {
        while (t < t_next)	/* evolve from t to t_next */
        {
            gsl_odeiv_evolve_apply (evolve_ptr, control_ptr, step_ptr,
                                    &my_system, &t, t_next, &h, y);
        }
        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); /* print at t=t_next */
    }

    /* all done; free up the gsl_odeiv stuff */
    gsl_odeiv_evolve_free (evolve_ptr);
    gsl_odeiv_control_free (control_ptr);
    gsl_odeiv_step_free (step_ptr);

    return 0;
}

/*************************** rhs ****************************/
/* 
   Define the array of right-hand-side functions y[i] to be integrated.
   The equations are:
   x' = v                  ==>  dy[0]/dt = f[0] = y[1]
   v' = -x + \mu v (1-x^2) ==>  dy[1]/dt = f[1] = -y[0] + mu*y[1]*(1-y[0]*y[0])
  
   * params is a void pointer that is used in many GSL routines
   to pass parameters to a function
*/
int
rhs (double t, const double y[], double f[], void *params_ptr)
{
    /* get parameter(s) from params_ptr; here, just a double */
    double mu = *(double *) params_ptr;

    /* evaluate the right-hand-side functions at t */
    f[0] = y[1];
    f[1] = -y[0] + mu * y[1] * (1. - y[0] * y[0]);

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}

/*************************** Jacobian ****************************/
/*
  Define the Jacobian matrix using GSL matrix routines.
  (see the GSL manual under "Ordinary Differential Equations") 
  
  * params is a void pointer that is used in many GSL routines
  to pass parameters to a function
*/
int
jacobian (double t, const double y[], double *dfdy,
          double dfdt[], void *params_ptr)
{
    /* get parameter(s) from params_ptr; here, just a double */
    double mu = *(double *) params_ptr;

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);

    gsl_matrix *m_ptr = &dfdy_mat.matrix;	/* m_ptr points to the matrix */

    /* fill the Jacobian matrix as shown */
    gsl_matrix_set (m_ptr, 0, 0, 0.0);	/* df[0]/dy[0] = 0 */
    gsl_matrix_set (m_ptr, 0, 1, 1.0);	/* df[0]/dy[1] = 1 */
    gsl_matrix_set (m_ptr, 1, 0, -2.0 * mu * y[0] * y[1] - 1.0); /* df[1]/dy[0] */
    gsl_matrix_set (m_ptr, 1, 1, -mu * (y[0] * y[0] - 1.0));     /* df[1]/dy[1] */

    /* set explicit t dependence of f[i] */
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}

Here is another cool example: https://github.com/cirosantilli/cpp-cheat/blob/6b4896846dfdca69e8d3f38fd3ede14c517438ef/gsl/ode.md