What would be an equivalent of "MaxSteps" using the GSL's ODE solver?

Question

I want to reproduce an ODE solver created using Mathematica with GSL.

Here is the Mathematica code which uses NDSolve:

result[r_] := NDSolve[{
    s'[t] == theta - (mu*s[t]) - ((betaA1*IA1[t] + betaA2*IA2[t] + betaB1*IB1[t] + betaB2*IB2[t]) +
                                  (betaA1T*TA1[t] + betaA2T*TA2[t] + betaB1T*TB1[t] + betaB2T*TB2[t])) * s[t] - 
                                 ((gammaA1*IA1[t] + gammaA2*IA2[t] + gammaB1*IB1[t] + gammaB2*IB2[t]) + 
                                  (gammaA1T*TA1[t] + gammaA2T*TA2[t] + gammaB1T*TB1[t] + gammaB2T*TB2[t])),

//... Some other equations



s[0] = sinit,IA1[0] = IA1init,IA2[0] = IA2init,
IB1[0] = IB1init,IB2[0] = IB2init,TA1[0] = TA1init,
TA2[0] = TA2init,TB1[0] = TB1init,TB2[0] = TB2init},
{s,IA1,IA2,IB1,IB2,TA1,TA2,TB1,TB2},{t,0,tmax},
MaxSteps->100000, StartingStepSize->0.1, Method->{"ExplicitRungeKutta"}];

Trying to get the exact equivalent using GSL:

int run_simulation() {
    gsl_odeiv_evolve*  e = gsl_odeiv_evolve_alloc(nbins);
    gsl_odeiv_control* c = gsl_odeiv_control_y_new(1e-17, 0);
    gsl_odeiv_step*    s = gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, nbins);
    gsl_odeiv_system sys = {function, NULL, nbins, this };
    while (_t < _tmax) {  //convergence check here
        int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &_t, _tmax, &_h, y);
        if (status != GSL_SUCCESS) { return status; }
    }
    return 0;
}

Where nbins is the number of equations given to the solver and _h the current step size.

I don't provide the equations themselves here, but the only way I found to limit the number of steps (as done with MaxSteps->100000 under Mathematica), is to adapt the first argument of the gsl_odeiv_control_y_new control feature. Here 1e-17 gives me something around 140000 steps...

Does anyone know a way to force the GSL's ODE solver to use a given maximum number of steps? As you probably understood, it is important to me to have results that I can really compare between those two tools.

Thanks for the help.

Answer 1

MaxSteps in Mathematica is only important when RK (Runge Kutta) gets stuck, and consequently fail to proper evolve your system. It does not fix the number of steps you want to take or the accuracy you need. Of course, higher accuracy demands lower step size which will imply more steps in a fixed interval. But my point is, unless you have a weird system where RK gets stuck and fails (and you would clear see the Mathematica error message in this case) or you set maxsteps to be ridiculous small, MaxSteps won't help you to proper compare mathematica and GSL.

To make a proper comparison you need to setup the same accuracy demands and control function in both programs. In fact, you can setup an arbitrary control function in GSL, besides the standard options, trough the API gsl_odeiv2_control_alloc and gsl_odeiv2_control_hadjust functions. You also must check what is the exactly stopping condition used in your Mathematica code.

Another option is to use a non-adaptive fixed step RK in both programs (in gsl you can call evolve the system with fix steps by calling gsl_odeiv2_driver_apply_fixed_step).

Last thing. 1e-17 seems to be an insane relative accuracy demand. Remember that roundoff errors usually does not allow RK to reach this level of accuracy. Actually roundoff errors is one of the things that can make RK to get stuck and/or make Mathematica/GSL to disagree with each other!!!! You should set accuracy to be > 1e-10.