仿真在運行ode45時掛起

Question

我想爲藻類羣體建立一個模型。這是我迄今爲止的代碼（所有代碼都是在線實例編寫的）。當我運行Solve_algaepop時，它只是掛了很長時間。

任何想法爲什麼？有什麼明顯的事情我做錯了嗎？方程來自一篇研究論文。

這是Solve_algaepop.m。在針對r1和r2的等式中，P10和P20被假定爲在algaepop_model.m中定義的值 P1 = x(1)和P2 = x(2)。我不知道如何當我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);