使用 GSL 的 ODE 求解器的“MaxSteps”等价物是什么？

Question

我想重现使用 Mathematica 和 GSL 创建的 ODE 求解器。

这是使用 NDSolve 的 Mathematica 代码：

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 是当前步长。

我没有在这里提供方程式本身，但我发现限制步数的唯一方法（就像在 Mathematica 下使用 MaxSteps->100000 所做的那样）是调整 gsl_odeiv_control_y_new 控制功能的第一个参数。这里1e-17 给了我大约 140000 步...

有谁知道强制 GSL 的 ODE 求解器使用给定的最大步数的方法？正如您可能理解的那样，对我来说，能够真正比较这两种工具的结果很重要。

感谢您的帮助。

Answer 1

Mathematica 中的MaxSteps 仅在 RK (Runge Kutta) 卡住并因此无法正确发展您的系统时才重要。它不固定您想要采取的步骤数或您需要的准确性。当然，更高的精度需要更低的步长，这意味着在固定间隔内有更多的步长。但我的观点是，除非您有一个奇怪的系统，其中 RK 卡住并失败（并且您会清楚地看到在这种情况下的 Mathematica 错误消息）或者您将 maxsteps 设置为荒谬的小，否则 MaxSteps 不会帮助您正确比较mathematica和 GSL。

要进行适当的比较，您需要在两个程序中设置相同的精度要求和控制功能。事实上，您可以在 GSL 中设置任意控制功能，除了标准选项外，还可以通过 API gsl_odeiv2_control_alloc 和 gsl_odeiv2_control_hadjust 函数。您还必须检查 Mathematica 代码中使用的确切停止条件是什么。

另一种选择是在两个程序中使用非自适应固定步骤 RK（在 gsl 中，您可以通过调用 gsl_odeiv2_driver_apply_fixed_step 来使用固定步骤来进化系统）。

最后一件事。 1e-17 似乎是一个疯狂的相对精度要求。请记住，舍入误差通常不允许 RK 达到这种准确度。实际上，舍入错误是可能使 RK 卡住和/或使 Mathematica/GSL 彼此不同意的事情之一！！！！您应该将精度设置为 > 1e-10。