代码之家 › 专栏 › 技术社区 › rocksNwaves

ODE45收敛到正确的曲线形状,但有错误的解决方案

numerical-analysis numerical-methods matlab

rocksNwaves · 技术社区 · 6 年前

提前谢谢你的帮助。我不是在寻找解决问题的明确方法,而是要指出我可能存在的明显错误。

我一直致力于在Matlab中求解一个非线性的一阶常微分方程组。在本研究中,对系统进行了数值求解: http://web.math.ku.dk/~moller/e04/bio/ludwig78.pdf

我一直在跟踪 documentation for ode45 ,并有运行的代码。

我已经做了所有的工作来理解和从头开始重建模型。我为一个班级的项目介绍了定性部分。我现在所做的是通过使用runge-kutta(或任何有效的方法)在Matlab中求解系统,使这个项目更进一步。最后,我想深入研究数值分析背后的理论,找出所选方法收敛的原因。

以下是数值求解系统的绘图,我正在尝试重新创建:

我发现我可以创建一个形状大致相同的情节,但有几个问题:

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

所以我要找的是这些差异的原因。我已经检查过我的odes系统和参数值很多次了,我的眼睛都模糊了。也许我在概念上遗漏了什么?

代码:

% 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

0 回复 | 直到 6 年前

Jamie Parkinson 6 年前

我看不出你的代码有什么问题。但我认为在产生这个数字的过程中有一些微妙之处,这在论文中没有很好的解释。

1)标度S轴(标签上写着“相对”)。我相信它们已经按k嫒s缩放了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

对初始条件有很大的敏感性。

进一步的调整可以使你更接近原始的数字,但我不确定这是否只是一个小技巧:在第一个方程中,用k_b替换k_b*y(2),没有这个,在减少之前,虫体密度会变得太大。新的情节如下。