[Help-gsl] Math error handling in ODE RK4

[Javier Andrés Mena Zapata]

Hi,

I currently have a model and I'm using the GSL RK4 implementation to
solve it, however it is not working as spected and I don't know how to
fix the problem. I implemented the "same" (I think) model in
Mathematica but the results are VERY differente... in GSL the result
function something like a decreasing exponential function but in
Mathematica it gives me a very smooth curve with maximums and
minimums.

The function that defines the ODE system is using an interpolation
function (created by me, but it's very simple and I'm sure it works
fine), and I think the problem is there.

Can you give a brief idea how can I fix the error?.

Here is a part of the code (in spanish, but several words are very
similar to the english meaning).

Thanks for your help.

struct informacion_problema {
 t_f_interpolacion f_interpolacion;
 vector_ciclico Hs;
 double rho;
 double r;
 double tau;
};

/* This function defines the ODE system */
int func (double t, const double y[], double f[],
     void* params)
{
 informacion_problema* mi_info = static_cast<informacion_problema*>(params);
 double H       = mi_info->f_interpolacion(t,mi_info->Hs);
 double rho     = mi_info->rho;
 double r       = mi_info->r;
 double tau     = mi_info->tau;
 double gamma_h = gamma(H,tau);

 f[0] = -H*y[0] +     rho*y[1] + gamma_h*y[2];
 f[1] =  H*y[0] - (rho+r)*y[1];
 f[2] =                 r*y[1] - gamma_h*y[2];

 return GSL_SUCCESS;
}

int main (void) {
 const gsl_odeiv_step_type* T = gsl_odeiv_step_rk4;
 double hrk4 = 1e-6;
 gsl_odeiv_step* s    = gsl_odeiv_step_alloc (T, 3);
 gsl_odeiv_control* c = gsl_odeiv_control_y_new (hrk4, 0.0);
 gsl_odeiv_evolve* e  = gsl_odeiv_evolve_alloc (3);
 informacion_problema mi_info;
 int tInterp; // tipo de interpolaci'on: ver m'as abajo
 double tinicial, tfinal, hpaso, y[3];

 HERE IS THE WHO READS THE PARAMS IN mi_info.
 USUALLY hpaso = 0.2

 gsl_odeiv_system sys = {func, NULL, 3, &mi_info};  // 3 is dimension
 double t = tinicial;
 double delta_t = (tfinal - tinicial);
 int periodos = static_cast<int>(round(delta_t/hpaso));

 WriteUtil::writeInt(periodos);

 for (int i = 1; i <= periodos; i++){
   double ti = i * delta_t / static_cast<double>(periodos) + tinicial;

   WriteUtil::writeDouble(t);
   WriteUtil::writeDouble(y[0]);
   WriteUtil::writeDouble(y[1]);
   WriteUtil::writeDouble(y[2]);
                        
   while (t < ti) {
     gsl_odeiv_evolve_apply (e, c, s,
                             &sys,
                             &t, ti, &hrk4,
                             y);
   }    
 }

 FREEING THE RESOURCES

}

--
Javier Andrés Mena Zapata
University of  Valle
Cali - Colombia