0

同等のバージョンの Matlab スクリプトを C に実装しました。ode45 を実行するために、GNU 科学ライブラリを使用しました。しかし、ode45 はバージョンごとに異なる出力を生成します。私はしばらく働いていますが、問題を見つけることができません。gsl_odeiv2_driver_apply_fixed_step を使用して、Matlab と同じ手順を実行します。

Matlab スクリプト

function ExpoGrowthEqn
% code to solve the exponential growth equation 
% dN/dt = r*N, N(0)=N0

% parameter values
r=0.2;
N0=10;

% numerical parameters
step=0.25;
tspan=[0:step:10];

% EXACT solution is: N(t)=N0 * exp(r*t)
for i=1:length(tspan)
    N(i,1)=N0*exp(r*tspan(1,i));
end
plot(tspan,N,'ob') %plots EXACT solution
hold on

% solve ODE using ode45
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
    [t y] = ode45( @growth_eqn, tspan, N0, options, r);
    plot(t,y,'-r') %plot APPROXIMATE solution

    solutions=[t N y]
end

function dy = growth_eqn(t, y, r)
N=y(1)
dy=r*N;
end

Cコード

void fixed_step(void) {

    PARAM parameters;
    parameters.N = 11500000;
    parameters.beta = 0.1;
    parameters.gamma = 1/2;   

    double I0=500.0/parameters.N;
    double R0=9000000.0/parameters.N;
    double S0=1-I0-R0;
    double y[3] = {S0, I0, R0};

    gsl_odeiv2_system sys = {func, NULL, 3, &parameters};
    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-12, 1e-12, 0);

    printf ("=== Initial values ===\n");
    printf ("y[0]=%6.15lf y[1]=%6.15lf y[2]=%6.15lf\n", y[0], y[1], y[2]);
    int length = 53;
    const double step = 0.25;
    double ti; double t = 0.0, tant = t;
    int i;
    for (i = 0; i <= length*(1/step); i++)
    {
        printf ("%.2f(%3.d) -> t=%.2f %6.15lf %6.15lf %6.15lf\n", (float)tant, i, (float)t, y[0], y[1], y[2]);
        tant = t;
        int status = gsl_odeiv2_driver_apply_fixed_step (d, &t, step/8, 8, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error, return value=%d\n", status);
        }
    }
    gsl_odeiv2_driver_free (d);        
}

Matlab 出力

S                   I                   R
------------------------------------------------------------------
0.217347826086956   0.000043478260870   0.782608695652174
0.217347603416592   0.000038578485741   0.782613818097667
0.217347405784518   0.000034229664122   0.782618364551360
0.217347230418583   0.000030370797298   0.782622398784119
0.217347074884551   0.000026948321844   0.782625976793605
...
0.217345850226451   0.000000000007068   0.782654149766481
0.217345850226413   0.000000000006234   0.782654149767354
0.217345850226380   0.000000000005510   0.782654149768110
0.217345850226351   0.000000000004888   0.782654149768761
0.217345850226327   0.000000000004356   0.782654149769316
0.217345850226307   0.000000000003903   0.782654149769790
0.217345850226289   0.000000000003514   0.782654149770197
0.217345850226273   0.000000000003172   0.782654149770554
0.217345850226259   0.000000000002859   0.782654149770881
0.217345850226245   0.000000000002556   0.782654149771198

C 出力

=== Initial values ===
y[0]=0.217347826086956 y[1]=0.000043478260870 y[2]=0.782608695652174
0.00(   ) -> t=0.00 0.217347826086956 0.000043478260870 0.782608695652174
0.00(  1) -> t=0.25 0.217347589196436 0.000043715151390 0.782608695652174
0.25(  2) -> t=0.50 0.217347351015482 0.000043953332344 0.782608695652174
0.50(  3) -> t=0.75 0.217347111537070 0.000044192810757 0.782608695652174
0.75(  4) -> t=1.00 0.217346870754133 0.000044433593694 0.782608695652174
...
51.00(205) -> t=51.25 0.217258883883215 0.000132420464611 0.782608695652174
51.25(206) -> t=51.50 0.217258162689554 0.000133141658272 0.782608695652174
51.50(207) -> t=51.75 0.217257437570519 0.000133866777307 0.782608695652174
51.75(208) -> t=52.00 0.217256708504771 0.000134595843055 0.782608695652174
52.00(209) -> t=52.25 0.217255975470856 0.000135328876970 0.782608695652174
52.25(210) -> t=52.50 0.217255238447202 0.000136065900623 0.782608695652174
52.50(211) -> t=52.75 0.217254497412123 0.000136806935703 0.782608695652174
52.75(212) -> t=53.00 0.217253752343811 0.000137552004014 0.782608695652174

2 番目の列の最後に、データに大きな違いがあります。Matlab ではほぼ 0 ですが、C コードではそうではありません。

Matlab タイムスタンプ

t =
                0
0.250000000000000
0.500000000000000
0.750000000000000
1.000000000000000
1.250000000000000
1.500000000000000
1.750000000000000
...
52.250000000000000
52.500000000000000
52.750000000000000
53.000000000000000
4

2 に答える 2

0

申し訳ありませんが、問題の原因となっている小さなコードが含まれていませんでした。不足している 4 行を追加するために、C コードが編集されています。問題が見つかりました:

// ガンマ = 0.0 parameters.gamma = 1/2;

結果として 0 を返します。明らかにこれは間違っています。コンマで数値を宣言すると解決します。

// ガンマ = 0.5 parameters.gamma = 1.0/2.0;

于 2013-04-05T08:12:53.050 に答える
0

指定した場合、Matlab は内部計算で固定ステップ サイズを使用しませんtspan(こちらを参照)。選択した時間に出力します。したがって、2 つの異なる計算を行っています。最初の GSL の例を読んで、固定出力時間で適応時間ステッピングを使用する方法を確認してください。

于 2013-04-04T15:30:45.917 に答える