如何在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着眼于在不同的初始条件下，系统的状态（这里你的三个不同的种群）如何演化。它不会告诉你任何有关这些解决方案的时间依赖性，只是它们相对于彼此的进展情况。