Что было бы эквивалентно "MaxSteps" с использованием решателя ODE GSL?

Question

Я хочу воспроизвести решатель 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 использовать заданное максимальное количество шагов? Как вы, наверное, поняли, для меня важно иметь результаты, которые я действительно могу сравнить между этими двумя инструментами.

Спасибо за помощь.

Answer 1

MaxSteps в Mathematica важна только тогда, когда RK (Runge Kutta) застревает и, следовательно, не может должным образом развить вашу систему. Он не определяет количество шагов, которые вы хотите сделать, или точность, которая вам нужна. Конечно, более высокая точность требует меньшего размера шага, что подразумевает большее количество шагов за фиксированный интервал. Но я хочу сказать, что если у вас нет странной системы, в которой RK застревает и дает сбой (и в этом случае вы бы увидели сообщение об ошибке Mathematica), или если вы установите для maxsteps нелепые малости, MaxSteps не поможет вам правильно сравнить mathematica и GSL.

Для правильного сравнения необходимо установить одинаковые требования к точности и функции управления в обеих программах. Фактически, вы можете установить произвольную функцию управления в GSL, кроме стандартных опций, через API gsl_odeiv2_control_alloc а также gsl_odeiv2_control_hadjust функции. Вы также должны проверить, какое именно условие остановки используется в вашем коде Mathematica.

Другим вариантом является использование неадаптивного фиксированного шага RK в обеих программах (в gsl вы можете вызвать evolve систему с фиксированными шагами, вызвав gsl_odeiv2_driver_apply_fixed_step).

Последнее. 1e-17 кажется безумным требованием относительной точности. Помните, что ошибки округления обычно не позволяют РК достичь этого уровня точности. На самом деле ошибки округления - это одна из вещей, которая может заставить RK застрять и / или заставить Mathematica/GSL не согласиться друг с другом!!!! Вы должны установить точность, чтобы быть> 1e-10.