如何在MATLAB中绘制非线性微分方程组的解? [英] How to plot solutions of system of nonlinear differential equations in MATLAB?

问题描述

我正在尝试遵循此研究论文.我正在尝试复制第20页上的图7中的解决方案图.我有图7的屏幕快照:

I am trying to follow this research paper. I am trying to duplicate the graph of the solutions located in Figure 7 on page 20. I have a screenshot of Figure 7:

我首先要重新创建左图.有问题的系统是我的dX.这是我在m文件中的文件:

I would first like to recreate the left picture. The system in question is what I have as dX. Here is what I have in an m-file:

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所拥有的.我在命令窗口中键入以下内容:

ode45 has the syntax: [T,Y] = solver(odefun,tspan,y0). I get the tspan from the picture I posted. My initial conditions are:S0 = 10^4; Rl0 = 0; Rh0 = 0, so this is what I have as y0. I type the following into the command window:

>>[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]);

过去30分钟，MATLAB一直很忙，我的笔记本电脑开始变得非常热.因此，在完成绘制之前，我无法进行绘图，而且我也不知道代码中是否存在任何错误.我想知道是否有更好的方法可以获取系统dX的解决方案.

MATLAB has been busy for the past 30 minutes and my laptop is starting to get really hot. So I can't get to plotting until it's finished and I don't know if there are any errors in my code. I was wondering if there was perhaps a better way that I could get the solutions of the system dX.

推荐答案

我对照纸检查了ODE，发现一个错误.根据论文第9页的最后一段，N = S + Rl + Rh.正确的代码:

I checked the ODE against the paper and found one mistake. According to the last paragraph of page 9 in the paper, N = S + Rl + Rh. A corrected code:

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).他们对相同的事物进行评估，但前者看起来更干净.

Note a couple of syntactic changes. First, since you do not actually use the time value in your ODE, you can leave a ~ in place of the t that you had in your function definition as a stand-in for the unused input. Second, you can use the notation 8e-3 instead of 8*10^(-3). They evaluate to the same thing, but the former looks a little cleaner.

未经处理的图，即T = 0，显示如下.

The plot for no treatment, i.e. T = 0, is shown below.

我首先尝试使用MATLAB中的一些刚性ODE求解器来求解方程，从而解决了ODE中的问题.有关更多详细信息，请参见 MATLAB ODE求解器文档.基本上，我用ode15s和ode23s对其进行了求解，发现该解决方案不稳定(人口激增到无穷大).如果您的解决方案未解决，这些其他求解器也是牢记的好工具.有时其他人中的一个会起作用，它会给您您想要的东西，或者表明您在其他地方遇到另一个问题.

I figured out the issue in the ODE by first trying to solve your equation with some of the stiff ODE solvers in MATLAB. See the MATLAB ODE solver documentation for more details. Basically I solved it with ode15s and ode23s and found that the solution was unstable (population went off to infinity). These other solvers are good tools to keep in mind if your solution is hanging; sometimes one of the others will work and it will either give you what you want, or indicate that you have another problem somewhere else.

注意:我可以肯定的是，您在此处使用的不是相像".您只是在寻找具有给定初始条件的ODE解决方案. 相像考察了给定的系统状态(此处是您的三个不同种群)不同的初始条件集.它没有显示任何有关这些解决方案的时间依赖性的信息，而没有显示它们相对于另一解决方案的演化方式.

Note: I am fairly certain that you are not using "phase portrait" correctly here. You are just looking for the solution of the ODE with respect to time with a given set of initial conditions. A phase portrait looks at how the states of the system (your three different populations here) evolve given different sets of initial conditions. It does not show you anything about the time-dependence of these solutions, just how they evolve relative to one-another.

这篇关于如何在MATLAB中绘制非线性微分方程组的解?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！