1

我想重现使用 GSL 创建的 ODEMathematica求解器。

这是使用的 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 步...

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

谢谢您的帮助。

4

1 回答 1

2

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

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

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

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

于 2013-10-07T20:37:58.163 回答