如何在MATLAB中繪製非線性微分方程組的解？

Question

我想關注此research paper。我試圖複製第20頁圖7中的解決方案圖。我有一個圖7的屏幕截圖：

我首先想重新創建左側圖片。有問題的系統是我的dX。以下是我在一個m文件：

function dX = CompetitionModel(t,X)  
    bs = 8*10^(-3); 
    bl = 4*10^(-3); 
    bh = 6.4*10^(-3); 
    N = bs + bl + bh; 
    K = 10^8; 
    m1 = 2*10^(-5); 
    m2 = 9*10^(-9); 
    p = 5*10^(-13); 
    I = 10^(-3); 
    T = 0; 
    a = 0; 
    dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3)); 
     X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3)); 
     X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))]; 
end 

ode45的語法是：[T,Y] = solver(odefun,tspan,y0)。我從我發佈的圖片中獲得了tspan。我的初始條件是：S0 = 10^4; Rl0 = 0; Rh0 = 0，所以這就是我的y0。我鍵入以下命令窗口：

>>[t,X1] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 
>>[t,X2] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 
>>[t,X3] = ode45('CompetitionModel', [0,45000], [10^4, 0, 0]); 

MATLAB一直忙於在過去的30分鐘，我的筆記本電腦開始變得非常熱。所以我不能在繪製完成之前進行繪圖，我不知道我的代碼中是否有任何錯誤。我想知道是否有更好的方式可以獲得系統dX的解決方案。

Answer 1

我檢查了ODE對紙張，發現一個錯誤。根據論文第9頁的最後一段，N = S + Rl + Rh。已更正的代碼：

function dX = CompetitionModel(~,X)  
    bs = 8e-3; 
    bl = 4e-3; 
    bh = 6.4e-3; 
    N = sum(X); % this line was incorrect 
    K = 1e8; 
    m1 = 2e-5; 
    m2 = 9e-9; 
    p = 5e-13; 
    I = 1e-3; 
    T = 0; 
    a = 0; 
    dX = [X(1) * (bs * (1 - N/K) - I - T - m1) - p * X(1) * (X(2) + X(3)); 
     X(2) * (bl * (1 - N/K) - I - a*T - m2) + m1 * X(1) + p * X(2) * (X(1) - X(3)); 
     X(3) * (bh * (1 - N/K) - I - a*T) + m2 * X(2) + p * X(3) * (X(1) + X(2))]; 

請注意一些語法更改。首先，由於實際上並未在ODE中使用時間值，因此您可以將~代替您在函數定義中使用的t作爲未使用輸入的替代品。其次，您可以使用符號8e-3而不是8*10^(-3)。他們評估的是同樣的事情，但前者看起來有點乾淨。

下面顯示了不處理的情節，即T = 0。

我首先試圖與一些在MATLAB僵硬ODE求解器來解決你的公式計算出在ODE問題。有關更多詳細信息，請參見MATLAB ODE solver documentation。基本上我用ode15s和ode23s解決了這個問題，發現解決方案不穩定（人口無限）。如果您的解決方案掛起，這些其他求解器是記住的好工具;有時候其中一個會工作，它會給你你想要的，或者表明你在別的地方還有其他問題。

注：我相當確定您在這裏沒有正確使用「相位肖像」。您只是在一組給定的初始條件下尋找ODE關於時間的解決方案。 A phase portrait着眼於在不同的初始條件下，系統的狀態（這裏你的三個不同的種羣）如何演化。它不會告訴你任何有關這些解決方案的時間依賴性，只是它們相對於彼此的進展情況。