1

GSLMathematicaを使用して作成した ODE ソルバーを再現したいと考えています。

を使用する 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現在のステップ サイズです。

ここでは方程式自体は提供しませんが、(Mathematica で行われたように)ステップ数を制限する唯一の方法は、制御機能MaxSteps->100000の最初の引数を適合させることです。gsl_odeiv_control_y_newここで1e-17は、約140000ステップが得られます...

GSL の ODE ソルバーに特定の最大ステップ数を使用させる方法を知っている人はいますか? おそらくご存じのとおり、これら 2 つのツールを実際に比較できる結果を得ることは、私にとって重要です。

助けてくれてありがとう。

4

1 に答える 1

2

MaxStepsIn Mathematica は、RK (Runge Kutta) が行き詰まってシステムを適切に進化させることができない場合にのみ重要です。実行したいステップ数や必要な精度を修正するものではありません。もちろん、精度が高いほど、ステップ サイズを小さくする必要があります。これは、一定の間隔でより多くのステップを実行することを意味します。しかし、私が言いたいのは、RK が動かなくなって失敗する奇妙なシステムを持っている場合 (この場合、Mathematica のエラー メッセージが表示されることは明らかです)、または maxsteps を途方もなく小さく設定しない限り、MaxSteps は mathematica を適切に比較するのに役立ちません。そしてGSL。

適切な比較を行うには、両方のプログラムで同じ精度要求と制御機能をセットアップする必要があります。gsl_odeiv2_control_alloc実際、APIと関数を介して、標準オプションに加えて、GSL で任意の制御機能をセットアップできますgsl_odeiv2_control_hadjust。また、Mathematica コードで使用されている正確な停止条件を確認する必要があります。

もう 1 つのオプションは、両方のプログラムで非適応固定ステップ RK を使用することです (gsl では、を呼び出して固定ステップでシステムを進化させることができますgsl_odeiv2_driver_apply_fixed_step)。

最後のこと。1e-17 は非常識な相対精度要求のようです。通常、丸め誤差によって RK がこのレベルの精度に達することはできないことに注意してください。実際、丸め誤差は、RK が動かなくなったり、Mathematica と GSL が一致しなくなったりする原因の 1 つです!!!! 精度を > 1e-10 に設定する必要があります。

于 2013-10-07T20:37:58.163 に答える