Я хочу воспроизвести решатель ODE, созданный с помощью Mathematica, с GSL< /а>.
Вот код Mathematica, в котором используется 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"}];
Попытка получить точный эквивалент с помощью 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;
}
Где nbins — количество уравнений, заданных решателю, а _h — текущий размер шага.
Я не привожу здесь сами уравнения, но единственный найденный мной способ ограничить количество шагов (как это делается с MaxSteps->100000 в системе Mathematica) — это адаптировать первый аргумент функции управления gsl_odeiv_control_y_new. Здесь 1e-17 дает мне что-то около 140000 шагов...
Кто-нибудь знает способ заставить решатель ODE GSL использовать заданное максимальное количество шагов? Как вы, наверное, поняли, для меня важно иметь результаты, которые я действительно могу сравнить между этими двумя инструментами.
Спасибо за помощь.
