GSL adaptive integration giving Invalid Pointer error

Question

I am trying to use the built-in RK2imp method for solving a slightly modified version of the van der pol problem given in the GSL sample examples for using adaptive solvers for integration. However, I keep getting this "error: invalid pointer". I am attaching my code here

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

#define ABTOL 1e-18

#define A 2.00861986087484313650940188
#define mu 1.0
#define TEND 6.6632868593231301896996820305

double eps_abs=ABTOL, eps_rel=0.0, tend=TEND;

int func (double t, const double y[], double f[], void *params)
{
    (void)(t); /* avoid unused parameter warning */
    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)
{
    (void)(t); /* avoid unused parameter warning */
    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 gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk2imp;

    gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, 2);
    gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
    gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (2);

    gsl_odeiv2_system sys = {func, jac, 2, NULL};

    double t = 0.0, t1 = TEND;
    double h = 1e-6;
    double y[2] = { A, 0.0 };

    while (t < t1)
    {
        int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &h, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error: %s\n", gsl_strerror (status));
            break;
        }

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

    gsl_odeiv2_evolve_free (e);
    gsl_odeiv2_control_free (c);
    gsl_odeiv2_step_free (s);
    return 0;
}
 

This is the output I get:

 error: invalid pointer

I am aware that I could be making a very simple (and silly) mistake but for some reason I am not able to see it. I would really appreciate any help! Thank you.

Answer 1

As I understand it, the inappropriate stepping function for the numerical integration method (Runge-Kutta). It’s better to use explicit methods here, you can even use gsl_odeiv2_step_rk2 stepping function, but the solution will converge for a very long time. OR you can use implicit integration methods, but you need to use step driver:

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

#define ABTOL 1e-18

#define A 2.00861986087484313650940188
#define mu 1.0
#define TEND 6.6632868593231301896996820305

double eps_abs=ABTOL, eps_rel=0.0, tend=TEND;

int func (double t, const double y[], double f[], void *params)
{
    (void)(t); /* avoid unused parameter warning */
    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)
{
    (void)(t); /* avoid unused parameter warning */
    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 gsl_odeiv2_step_type * T = gsl_odeiv2_step_msadams;
    // with "gsl_odeiv2_step_rk2imp" stepping function the solution
    // converges VERY SLOWLY and may be "failure" error at the end

    gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, 2);
    gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
    gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (2);

    gsl_odeiv2_system sys = {func, jac, 2, NULL};

    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(&sys, T, 1e-6, 1e-6, 1e-6 );
    gsl_odeiv2_step_set_driver(s, d);

    double t = 0.0, t1 = TEND;
    double h = 1e-6;
    double y[2] = { A, 0.0 };

    while (t < t1)
    {
        int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &h, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error: %s\n", gsl_strerror (status));
            break;
        }

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

    gsl_odeiv2_evolve_free (e);
    gsl_odeiv2_control_free (c);
    gsl_odeiv2_step_free (s);
    return 0;
}