什麼是使用GSL的ODE求解器的「MaxSteps」的等價物？

我想重現使用Mathematica與GSL創建的ODE求解器。什麼是使用GSL的ODE求解器的「MaxSteps」的等價物？

下面是使用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解算器使用給定的最大步數？正如你可能理解的那樣，對我而言，能夠真正比較這兩種工具的結果對我來說很重要。

感謝您的幫助。

來源

2013-10-06 Gauthier Boaglio

MaxSteps在Mathematica中只有在RK（Runge Kutta）卡住時才重要，因此無法正確演進您的系統。它不會修復您想要採取的步驟數量或您需要的準確度。當然，更高的精確度要求更小的步長，這意味着在固定的時間間隔內會有更多的步驟。但是我的觀點是，除非你有一個奇怪的系統，RK卡住並失敗（在這種情況下你會清楚看到Mathematica錯誤信息），或者你設置的maxsteps很小，MaxSteps不會幫助你正確比較數學和GSL。

爲了進行適當的比較，您需要在兩個程序中設置相同的精度要求和控制功能。實際上，除了標準選項，您還可以通過API gsl_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-07 20:37:58

非常感謝。在進一步挖掘之後，我對你提到的有關'MaxSteps'的相同結論......關於你的其他觀點，這聽起來很有趣，我會嘗試這些可能性（我不知道'gsl_odeiv2_control_hadjust'，我們也不能使用GSL中的固定步驟）。那謝謝啦！還有一個問題：函數來自哪裏？_ odeiv'2' _...函數來自哪裏？我沒有他們。我的GSL版本是否太舊了？ –

yeap ..你需要更新gsl到1.15或1.16來獲得odeiv2 –

好的，再次感謝Vinicius。乾杯;） –

什麼是使用GSL的ODE求解器的「MaxSteps」的等價物？

回答

相關問題