0

我在 C 中实现了一个 Matlab 脚本的等效版本。为了运行 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

在第二列的末尾，数据有很大不同。在 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

tspan当您指定（阅读此处）时，Matlab 在其内部计算中不使用固定步长。它在那些选定的时间输出。因此，您正在进行两种不同的计算。阅读第一个 GSL 示例，了解如何使用具有固定输出时间的自适应时间步进。

于 2013-04-04T15:30:45.917 回答
0

抱歉，我没有包含导致问题的一小段代码。C 代码已被编辑以添加缺少的 4 行。发现问题：

// gamma = 0.0 参数.gamma = 1/2;

结果返回 0。显然这是错误的。用逗号声明数字可以解决它。

// gamma = 0.5 参数.gamma = 1.0/2.0;

于 2013-04-05T08:12:53.050 回答