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];
}
}_______________________________________________ Help-gsl mailing list [EMAIL PROTECTED] http://lists.gnu.org/mailman/listinfo/help-gsl
