[Help-gsl] ODE solver gets stuck

[Paul Schneider]

I'm trying to solve a system of ODE's with GSL. I have to say that I 
have no experience with ODE's and especially not with numerically 
nolving them. Therefore please don't hold back even with very basic hints.

My problem: a 4 dimensional ODE can be solved rapidly, only for some 
parameter constellations GSL seems to get stuck in some loop and I have 
to manually kill the program. I looked at the output during such an 
incident and it seems that at some point t from below doesn't increase 
anymore, but the function values grow infinitely. I have tried all the 
solvers, with and without Jacobian. The following code is what I think 
is the relevant portion of my C++ program using GSL. TinyVector (PVec) 
is from the Blitz library.

The function odeSolve below is called thousands of times for different 
(global) parameters KQd_ThQd_, KQd_. I am very grateful to all hints how 
I can get this to work. The ODE itself has not been reported to be 
especially difficult to handle.

Thanks,

Paul

int odeQd(double t, const double y[], double f[],    void *params)
{
   PVec result(0);
   PVec result4(0);
   double result2 = 0;
   double result3 = 0;
   double check;
   int status = GSL_SUCCESS;
   for ( int i = 0; i < M::N; ++i){
       result2 += KQd_ThQd_[i] * y[i+1]; // first diff first expression
       result3 += y[i+1] * y[i+1] * alpha_[i]; // first diff second 
expression
       // this is a vector expression
       result4 +=  (y[i+1] * y[i+1]) * beta_[i]; // second vector ODE 
second expression
       for ( int j = 0; j < M::N; ++j){
           result[j] -= y[i+1] * KQd_(i, j); // second vector ODE first 
expression
       }
   }
   // Handle errors
   check = (result2) + 0.5 * result3 - delD0_;
   if (!gsl_finite(check)){
       GSL_ERROR ("NaN or infinity encountered in coeff a", GSL_ERANGE);
       status = GSL_ERANGE;
   }
   f[0] = check;
   result += (0.5 * result4 - delD_);
   for ( int i = 0; i < M::N; ++i ){
       check = result[i];
       if (!gsl_finite(check)){
           GSL_ERROR ("NaN or infinity encountered in coeff B", 
GSL_ERANGE);
           status = GSL_ERANGE;
       }
       f[i+1] = check;
   }

   return status;
}

const gsl_odeiv_step_type * ODE_T_   = gsl_odeiv_step_rkf45;
gsl_odeiv_step * ODE_s_;
gsl_odeiv_control * ODE_c_;
gsl_odeiv_evolve * ODE_e_;

gsl_odeiv_system domSys = {odeQd, jacQd, M::N+ 1, 0};



double const ODE_absErr_ = 0.0;
double const ODE_relErr_ = 1e-6;

   ODE_s_ = gsl_odeiv_step_alloc (ODE_T_, M::N+1);
   ODE_c_ = gsl_odeiv_control_y_new (1e-6, 0.0);
   ODE_e_ = gsl_odeiv_evolve_alloc (M::N+1);

void odeSolve(double tau, double& a, PVec& vec, IDN ide)
{
   gsl_odeiv_step_reset(ODE_s_);
   gsl_odeiv_evolve_reset(ODE_e_);

   blitz::TinyVector<double, M::N+1> initial(0);
   double t = 0.0;
   double h = 1e-6;
   while (t < tau)
   {
       int status = gsl_odeiv_evolve_apply (ODE_e_, ODE_c_, ODE_s_,
                                            &sys,
                                            &t, tau,
                                            &h, initial.dataFirst());

       if (status != GSL_SUCCESS){
           break;
       }
   }
   a = initial[0];
   for( int i = 1; i <= M::N ; ++i ){
       vec[i-1] = initial[i];
   }
}