ode45 收敛到正确的曲线形状，但解错误

Question

在此先感谢您的帮助。我不是在为我的问题寻找明确的解决方案，而是要指出我可能明显的错误。

我一直致力于在 MATLAB 中求解非线性一阶 ODE 系统。该系统在本研究中被数值求解：http://web.math.ku.dk/~moller/e04/bio/ludwig78.pdf

我一直在关注 documentation for ode45，并且有运行的代码。

我已经完成了从头开始理解和重新创建模型的所有工作。我介绍了 class 项目的定性部分。我现在正在做的是通过使用 runge-kutta（或任何有效的方法）在 MATLAB 中求解系统，使该项目更进一步。最后，我想深入探讨数值分析背后的理论，找出所选方法收敛的原因。

这是我正在尝试重新创建的数值求解系统的图：

我发现我可以创建一个形状大致相同的图，但是有几个问题：

发生变化的时间尺度是上图的三倍。
函数值的范围大错特错。
仅当我将初始条件调整为与上面 t=0 附近显示的内容明显不同。

所以我正在寻找的是这些差异的原因。我已经多次检查我的 ODE 系统和参数值，以至于我的眼睛都模糊了。也许我在概念上遗漏了什么？

代码：

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

tspan = [0 200];
init = [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
subplot(3,1,1);
plot(t,Y(:,1),'b');
title('Budworm Density');

subplot(3,1,2)
plot(t,Y(:,2),'g');
title('Branch Density');

subplot(3,1,3);
plot(t,Y(:,3),'r');
title('Foliage Condition');

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*y(3)));
         r_e*y(3)*(1 - (y(3)/k_e)) - p*y(1)/y(2)
        ];

end

Answer 1

我没有发现您的代码有任何问题。但我认为生成该图有一些微妙之处，论文中没有很好地解释。

1) S轴被缩放（在标签中写着'relative'）。我相信他们已经将 S 缩小了 k_s。我认为您还需要缩放参数 p（设置 p = p*k_s），否则 E 等式中的最后一项将很小，并且 E 人口不会在所需的时间尺度内减少。

2) 我认为他们一定对 E 施加了一些下限，以避免被 0 除。你可以在图中看到 E->0 首先，但是在你的 S 方程中，如果发生这种情况，那么你将除以 0，求解器不会收敛。

将这些放在一起，对您的代码进行以下轻微修改会产生与论文中更相似的结果：

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

% Scale p with k_s
p = p*k_s;

tspan = [0 50]; % [0 200];
init = [1e-16 0.075*k_s 1]; % [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);

% To scale before plotting, so everything fits on a 0->1 y axis.
maxB = 500;
S_scale = k_s;

figure('Position', [200 200 1000 600]);
hold on;
plot(t,Y(:,1)/maxB,'b');
plot(t,Y(:,2)/(S_scale),'g');
plot(t,Y(:,3),'r');

ylim([0, 1]);

hold off;

box on;

legend({['Budworm Density, B / ', num2str(maxB)], 'Branch Density, S / 0.75', 'Foliage Condition, E'}, ...
    'Location', 'eastoutside')

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)

% Place lower limit on E
E = max(y(3), 1e-5);

dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*E));
         r_e*E*(1 - (E/k_e)) - p*y(1)/y(2)
        ];

end

对初始条件非常敏感。

进一步调整可以使您更接近原始图形，但我不确定这是否只是一个 hack：在第一个等式中，将 k_b*y(2) 替换为 k_b。没有这个，Budworm 密度在下降之前变得太大。新剧情如下。

ode45 收敛到正确的曲线形状，但解错误

ode45 converges to correct curve shape, but with wrong solution

matlab

numerical-analysis

numerical-methods