GSL ODE solver returns -nan although same ODE with same parameters is being solved in python

Question

I use python to solve ODEs using scipy.integrate.odeint. Currently, I am working on a small project where I am using gsl in C++ to solve ODEs. I am trying to solve an ODE but the solver is returning -nan for each time point. Following is my code:

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

struct param_type {
    double k;
    double n;
    double m;
    double s;
};


int func (double t, const double y[], double f[], void *params)
{
  (void)(t); /* avoid unused parameter warning */

  struct param_type *my_params_pointer = (param_type *)params;
  double k  = my_params_pointer->k;
  double n  = my_params_pointer->n;
  double m  = my_params_pointer->m;
  double s  = my_params_pointer->s;

  f[0] = m*k*pow(s,n)*pow((y[0]/(k*pow(s,n))),(m-1)/m);

  return GSL_SUCCESS;
}


int * jac;


int main ()
{
  struct param_type mu = {1e-07, 1.5, 0.3, 250};
  gsl_odeiv2_system sys = {func, NULL, 1, &mu};

  gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
  int i;
  double t = 0.0, t1 = 10.0;
  double step_size = 0.1;
  double y[1] = { 1e-06 };

  gsl_vector *time = gsl_vector_alloc ((t1 / step_size) + 1);
  gsl_vector *fun_val = gsl_vector_alloc ((t1 / step_size) + 1);

  for (i = 1; i <= t1/step_size; i++)
    {
      double ti = i * t1 / (t1 / step_size);
      int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

      if (status != GSL_SUCCESS)
        {
          printf ("error, return value=%d\n", status);
          break;
        }

      printf ("%.5e %.5e\n", t, y[0]);
      gsl_vector_set (time, i, t);
      gsl_vector_set (fun_val, i, y[0]);
    }

  gsl_vector_add(time, fun_val);

  {
     FILE * f = fopen ("test.dat", "w");
     gsl_vector_fprintf (f, time, "%.5g");
     fclose (f);
  }

  gsl_odeiv2_driver_free (d);
  gsl_vector_free (time);
  gsl_vector_free (fun_val);

  return 0;
}

As mentioned here, I don't need jacobian for an explicit solver that's why I passed NULL pointer for the jac function. When I run the above code, I get -nan values at all time points.

To cross-check, I wrote the program in python which has the same function and same parameters, solved using scipy.integrate.odeint. It runs and provides a plausible answer. Following my python code:

import numpy as np
from scipy.integrate import odeint

def nb(y, t, *args):
    k = args[0]
    n = args[1]
    m = args[2]
    s = args[3]
    return m*k*s**n*(y/(k*s**n))**((m-1)/m)

t = np.linspace(0,10,int(10/0.1))
y0 = 1e-06

k = 1e-07
n = 1.5
m = 0.3
s = 250

res = odeint(nb, y0, t, args=(k,n,m,s)).flatten()
print(res)

Could anyone please help me figure out, what I am doing wrong in the C++ code using GSL for solving the ODE?

Answer 1

Your problem is here:

f[0] = m*k*pow(s,n)*pow((y[0]/(k*pow(s,n))),(m-1)/m);

As the solver proceeds, it may want to sample negative values of y[0]. In Python this makes no problem, in C++ it produces NANs.

To handle this, you can mimic Python's behavior:

  auto sign = (y[0] < 0) ? -1.0 : 1.0;
  f[0] = sign*m*k*pow(s,n)*pow((std::abs(y[0])/(k*pow(s,n))),(m-1)/m);

or even set sign effectively to 1:

  f[0] = m*k*pow(s,n)*pow((std::abs(y[0])/(k*pow(s,n))),(m-1)/m);

After all, raising negative values to noninteger powers is an error unless one considers complex numbers, which is not the case.

Please notice that y[0] was secured with std::abs.