matlab - 运行 ode45 时模拟挂起

Question

我想为藻类种群建立模型。这是我到目前为止的代码(全部来自在线示例)。当我运行 Solve_algaepop 时，它会挂起很长时间。

有什么想法吗？有什么明显的我做错了吗？这些方程式来自一篇研究论文。

这是 Solve_algaepop.m。在 r1 和 r2 的等式中，P10 和 P20 应该是值P1 = x(1) 和 P2 = x(2) 在 algaepop_model.m 中定义。当我在 Solve_algaepop.m

时，我不知道如何访问这些值

% Initial conditions
P10 = 560000; %from Chattopadhyay; estimated from graph
P20 = 250000; %same as above
Z0 = 280000; %
N0 = 0.6; %from Edwards

%some variables that the expressions of the parameters use
lambda = .6;
mu = .035;
k = 0.05;

%define parameters (start with estimates from Edwards paper):
r1 = (N0/(.03+N0))*((.2*P10)/(.2 + .4*P10)); 
r2 = (N0/(.03+N0))*((.2*P20)/(.2 + .4*P20));
a = Z0*((lambda*P10^2)/(mu^2 + P10^2));%G1: zooplankton growth function from Edwards paper
% m1 = .15; %r in Edwards paper
m1 = .075; % q in Edwards
m2 = .15;% r in Edwards paper
m3 = .15; % r in Edwards paper
d = 0.5;
cN = k;%*(N-N0);

par = [r1 r2 a m1 m2 m3 d cN]; % Creates vector of parameter values to pass to the ode solver

tspan = 0:1:300; %(Note: can also use the function linspace)

x0 = [P10 P20 Z0 N0]; % Creates vector of initial conditions

[t,x] = ode45(@algaepop_model,tspan,x0,[],par);
plot(t,x)

这里是 algaepop_model.m

function dxdt = algaepop_model(t,x,par)

P1 = x(1);
P2 = x(2);
Z = x(3);
N = x(4);

r1 = par(1);
r2 = par(2);
a = par(3);
m1 = par(4);
m2 = par(5);
m3 = par(6);
d = par(7);
cN = par(8);

dxdt = zeros(4,1);

dxdt(1) = r1*N*P1 - m3*P1 - a*P1*Z;
dxdt(2) = r2*N*P2 - a*P2*Z - m2*P2;
dxdt(3) = a*P2*Z + a*P1*Z - m1*Z;
dxdt(4) = d*m2*P2 + d*m1*Z + d*m3*P1 + cN - r2*N*P2 - r1*N*P1;

end

感谢您的帮助。

Answer 1

让我们调试。您可以做的最简单的事情之一是在您的集成函数 algaepop_model 中打印出 t 和 x。一旦执行此操作，您可能会注意到发生了什么:ode45 正在采取非常小的步骤。它们的顺序为 1.9e-9。使用如此小的步骤，模拟到 t = 300 将花费很长时间(如果您在每一步都打印出内容，甚至会更长)。

这可能是由于初始条件选择不当、缩放比例或维度化不当、导致方程式错误的拼写错误，或者仅仅是因为您对特定问题使用了不合适的求解器(和/或公差)。我无法真正解决前两种情况，必须假设您没有任何错误。因此，在这种情况下，您实际上拥有一个 stiff system。在这种情况下，ode45 并不是一个特别好的选择。只需将求解器更改为 ode15s几乎立即产生以下情节:

如您所见，情节的初始部分在短时间内发生了非常大的变化。如果放大，您会看到在第一个时间单位内发生了巨大的尖峰(您可能会输出更多时间点，或者只是让 tspan = [0 300])。一些状态变量变化很快，而另一些状态变量变化得比较缓慢。如此高的频率和时间尺度的差异是刚性系统的标志。我建议，除了确认您的代码正确之外，您还可以尝试通过 odeset 调整集成公差。 .确保更严格的公差产生质量相似的结果。您也可以尝试 other stiff solvers in the ODE suite如果你愿意。

最后，通过函数句柄本身而不是您自己的方式来传递参数更有效和最新。方法如下:

[t,x] = ode15s(@(t,x)algaepop_model(t,x,par),tspan,x0);