c - What would be an equivalent of "MaxSteps" using the GSL's ODE solver? -
i want reproduce ode solver created using mathematica
gsl.
here mathematica code 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])), //... 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 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, }; 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
number of equations given solver , _h
current step size.
i don't provide equations here, way found limit number of steps (as done maxsteps->100000
under mathematica), adapt first argument of gsl_odeiv_control_y_new
control feature. here 1e-17
gives me around 140000 steps...
does know way force gsl's ode solver use given maximum number of steps? understood, important me have results can compare between 2 tools.
thanks help.
maxsteps
in mathematica important when rk (runge kutta) gets stuck, , consequently fail proper evolve system. not fix number of steps want take or accuracy need. of course, higher accuracy demands lower step size imply more steps in fixed interval. point is, unless have weird system rk gets stuck , fails (and clear see mathematica error message in case) or set maxsteps ridiculous small, maxsteps won't proper compare mathematica , gsl.
to make proper comparison need setup same accuracy demands , control function in both programs. in fact, can setup arbitrary control function in gsl, besides standard options, trough api gsl_odeiv2_control_alloc
, gsl_odeiv2_control_hadjust
functions. must check stopping condition used in mathematica code.
another option use non-adaptive fixed step rk in both programs (in gsl can call evolve system fix steps calling gsl_odeiv2_driver_apply_fixed_step
).
last thing. 1e-17 seems insane relative accuracy demand. remember roundoff errors not allow rk reach level of accuracy. roundoff errors 1 of things can make rk stuck and/or make mathematica/gsl disagree each other!!!! should set accuracy > 1e-10.
Comments
Post a Comment