微分方程模块

引例

列出方程组

image-20240129132038919

速度公式

一个导弹头对准的等式（斜率相等）

image-20240129132523940
image-20240129132602724

求解方程组

image-20240129132622152

概述

image-20240129132920565
image-20240129133016467

Matlab求微分方程的解析解

求解析解

image-20240129133211310
image-20240129133240995 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%% 例1
clear;clc
dsolve('y-Dy=2*x','x')  % 这里要指定自变量为x
% 2*x + C1*exp(x) + 2  (这里的C1表示任意常数，有时候也会出现C2 C3等)
dsolve('y-Dy=2*x')  % 如果不指定自变量的话，会默认自变量为t，x会看成一个常数
% 2*x + C2*exp(t)

% 注意：最新版本的matlab会逐渐淘汰上面那种写法（虽然我个人觉得上面的写法更方便）
% 下面这种写法是新版的matlab推荐的方式（和我们上一讲符号运算中解方程的写法类似）
syms y(x)
eqn = (y - diff(y,x) == 2*x);    % 注意原来方程中的“=”一定要改成“==”
dsolve(eqn)

%% 如果微分方程中还有其他的未知参数怎么办？
% 方法1
dsolve('y-Dy=a*x','x')  % a是一个未知的参数
% 方法2
syms y(x) a
eqn = (y - diff(y,x) == a*x);  
dsolve(eqn)

%% 例2 
% 方法1
dsolve('y-Dy=2*x','y(0)=3','x')
% 2*x + exp(x) + 2
% 方法2
syms y(x)
eqn = (y - diff(y,x) == 2*x);  
cond = (y(0) == 3);
dsolve(eqn,cond)
% 2*x + exp(x) + 2


%% 例3
% 方法1
dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 
% 3*sin(5*x)*exp(-2*x)
% 方法2
syms y(x)
eqn = (diff(y,x,2) + 4 *diff(y,x) + 29*y  == 0);  
Dy = diff(y,x); % 定义变量Dy为y的一阶导数
cond = [(y(0) == 0) ,(Dy(0) ==15)] ; % 有两个条件，可以写到一个向量中保存
dsolve(eqn,cond)
% 3*sin(5*x)*exp(-2*x)

%% 例4
% 方法1
[x,y,z] = dsolve('Dx=2*x-3*y+3*z+t','Dy=4*x-5*y+3*z+t','Dz=4*x-4*y+2*z+t','t') 

% 方法2
syms x(t) y(t) z(t)
eqn1 = (diff(x,t)  == 2*x-3*y+3*z+t); 
eqn2 = (diff(y,t)  == 4*x-5*y+3*z+t); 
eqn3 = (diff(z,t)  == 4*x-4*y+2*z+t); 
eqns = [eqn1 eqn2 eqn3];
[x,y,z] = dsolve(eqns)
% x = exp(2*t)*(C2- (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t)
% y = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t) + C4*exp(-2*t)
% z = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C4*exp(-2*t)
mupad  % 最新版本matlab可能会报错，将计算结果复制到里面，使结果可读。
% 如果新版matlab用不了mupad的话，可以使用更新13中介绍到的实时脚本
simplify(y)  % simplify函数可以简化表达式
latex(y) % 转换成latex代码，复制到Axmath或者word自带的公式编辑器（低版本不知道支不支持）
% 如果太过于复杂的话可能会报错，大家可以自己测试

%% 不是所有的微分方程都可以，导弹追击那一题就没有解析解
% 假设 v=100
[x,y] = dsolve('Dx = 3*100*(20+sqrt(2)/2*100*t-x)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','Dy = 3*100*(sqrt(2)/2*100*t-y)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','x(0)=0,y(0)=0','t')  
% 警告: Explicit solution could not be found.

Matlab求解微分方程的数值解

image-20240129133711729

函数如下

image-20240129134021044

标准形式

image-20240129134123447

一般用ode45和ode15s

image-20240129134325396
image-20240129134427546

判断标准：刚性和非刚性的判断

image-20240129134449713

基本都是非刚性问题，因此一般用ode45，百分之90

例题一：

image-20240129134704414 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%% 调用ode45函数求解微分方程df1，自变量为x，范围为[0,2],  初始值y(0)=3 ; 因变量为y
clear;clc
[x,y] = ode45('df1',[0,2],3);  % [x,y] = ode45(@df1,[0,2],3);
% [x,y] = ode23('df1',[0,2],3);
% [x,y] = ode113('df1',[0,2],3);
% [x,y] = ode15s('df1',[0,2],3); % 刚性
figure(1)
plot(x,y,'r*-')  % 画出x和y的函数图像，用红色的直线和*标记

%%  下面我们直接画出微分方程的解析解的图像进行对比
hold on
x = 0:0.01:2;
y = exp(x)+2*x+2;
plot(x,y,'b-')  % 蓝色的直线
% 从图中可以看出，ode45函数得到的数值解的精度很高。

%%  设定相对误差和绝对误差，这样可以提高微分方程数值解的精度
options = odeset('reltol',1e-4,'abstol',1e-8);
[x,y] = ode45('df1',[0,2],3,options);

%%  如果觉得x的间隔不够小，我们可以指定要求解的位置
[x,y] = ode45('df1',[0:0.001:2],3,options);





function dy = df1(x,y)
    % 微分方程：y-y'=2x（函数名称可以任意取）
    dy = y - 2*x; % 写成标准形式 y' = y - 2x 
    % 注意函数的返回值一定是因变量y的一阶导数
    % 函数的输入有两个，分别是自变量x和因变量y
end

例题二：

image-20240129135125772 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%% 在真正比赛中，我们往往会倾向于先去看有无解析解，如果没有解析解再考虑数值解
% [y1 y2 y3] = dsolve('Dy1=y2*y3','Dy2=-y1*y3','Dy3=-0.51*y1*y2','y1(0)=0,y2(0)=1,y3(0)=1','x')

%% 调用ode45函数求解微分方程df2，自变量为x，范围为[0,4*pi] ; 因变量为y1,y2和y3，初始值: y1(0)=0,y2(0)=y3(0)=1 
[x, y] = ode45('df2', [0 4*pi], [0 1 1]);  % 这里的y是一个有3列的矩阵哦！
plot(x, y(:,1), 'o', x, y(:,2), '*', x, y(:,3), '+') 
legend('y1','y2','y3')  % 加上标注
axis([0, 4*pi, -inf, +inf])  % 设置横坐标范围为0-4pi，纵坐标范围不需要设置，写成-inf到+inf

function dy = df2(x,y) 
       % 注意哦，x是自变量，y是因变量，由y1,y2,y3组成 
       dy = zeros(3,1);  % 初始化用来储存因变量一阶导数的列向量（不能写成行向量哦）
       dy(1) = y(2) * y(3);
       dy(2) = -y(1) * y(3);
       dy(3) = -0.51 * y(1) * y(2);
%     上面四行可以写成一行：   dy = [ y(2) * y(3);   -y(1) * y(3);  -0.51 * y(1) * y(2)]
end

例题三：

image-20240129135929028
image-20240129140127574
image-20240129140153081 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% 我们先看这个微分方程有没有解析解
dsolve('(1+x*x)*D2y=2*x*Dy','y(-2)=3,Dy(-2)=4','x')
x = -2:0.01:2;
y = (4*x.*(x.^2 + 3))/15 + 101/15;
plot(x,y,'b-')
% 事实上有了解析解后我们就完全不需要数值解了，下面我们只是为了演示怎么进行求解~
hold on % 继续在上面作图
[x,y]=ode45('df3',[-2,2],[3,4]);  % 求出这个微分方程df3的数值解
plot(x,y(:,1),'r*') % 注意，y变量有两列，第一列是y1(我们要求的y)，第二列是y2(y的一阶导数) 


function dy=df3(x,y) 
       dy=zeros(2,1);  % 初始化用来储存因变量一阶导数的矩阵
       dy(1)=y(2);
       dy(2)=(2*x)/(1+x*x)*y(2);
end

例题4:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% 我们先看这个微分方程有没有解析解
dsolve('D2y=1000*(1-y^2)*Dy-y','y(0)=2,Dy(0)=0','x')  % 警告: Explicit solution could not be found. 

% 下面计算数值解
% 如果使用ode45函数会发现计算的非常慢，matlab一直显示正忙(windows电脑在命令行窗口按Ctrl+C可以终止运行)
% [x,y]=ode45('df4',[0,3000],[2,0]);  % 求出这个微分方程df4的数值解
% 所以我们可以使用刚性问题的函数ode15s对其求解
[x,y]=ode15s('df4',[0,3000],[2,0]);  % 求出这个微分方程df4的数值解
plot(x,y(:,1),'*') % 注意，y变量有两列，第一列是y1(我们要求的y)，第二列是y2(y的一阶导数)

function dy=df4(x,y)
    dy=zeros(2,1);
    dy(1)=y(2);
    dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
end


例题5：导弹追踪问题

image-20240129140449477 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
clear; clc
options = odeset('reltol',1e-4,'abstol',1e-8);  %  设定相对误差和绝对误差，这样可以提高微分方程数值解的精度
% 下面的[0,0.1]表示时间t的范围，因为我们导弹速度假设的是200，所以这里的范围给小点
[t,y]=ode45('df5' ,[0,0.1],[0 0], options); % y是因变量，第一列为导弹运行的横坐标，第二列为导弹运行的纵坐标
% [t,y]=ode45('df5' ,[0:0.0001:0.1],[0 0], options); % y是因变量，第一列为导弹运行的横坐标，第二列为导弹运行的纵坐标
x = y(:,1);  y =y(:,2);
plot(x, y,'*','MarkerSize',1)  % 画出导弹的运行轨迹
hold on
plot([20,30],[0,10]) % 画出B船的运行轨迹: x-y-20=0
hold on
% 下面我们找到导弹与B船相撞的点（由于Matlab浮点数计算的原因，距离足够近时即可认为相撞）
n =length(t);  % 找到Matlab计算微分方程的数值解时一共有多少个时间点
d = 0;  % 初始化导弹飞行的距离
for i = 1:n % 开始循环
    dd = abs(x(i) - y(i) - 20) / sqrt(2);  % 利用点到直线的距离公式计算导弹和船的距离
    if dd < 0.001 % 如果这个距离足够小了，我们就认为相撞了，但再此之前别忘了判断导弹是否达到了有效射程
        for k = 2:i
            d = sqrt((x(k)-x(k-1))^2+(y(k)-y(k-1))^2) + d;  % 以直代曲的思想求曲线的长度，即导弹飞行的距离
        end
        if d <= 50 % 导弹的有效射程为50个单位,如果没有达到50单位
            disp(['导弹飞行',num2str(d),'单位后击中B船'])
            disp(['导弹飞行的时间为',num2str(t(i)*60 ),'分钟']) % 输出导弹击中B船的时间(转换为分钟)
            disp('击中点的坐标：') ;  disp([x(i),y(i)])  % 输出导弹击中B船的坐标
            plot(x(i),y(i),'r*');
            text(x(i)+0.5,y(i)+0.1,'击中点')
        end
        break; % 跳出循环
    end
end
if d >50 || dd >= 0.001 % 如果射程大于50或导弹与B船的距离始终都没有小于0.001（这个数需要根据实际情况调整）
    disp('导弹没有击中B船');
end

t(i) * 200 * 3   % 更快计算导弹飞行距离的公式：速度*时间
% 得到的结果和我们上面以直代曲的结果很接近


% 下面是以前使用蒙特卡罗模拟得到的解，大家可以对比下：
% 导弹飞行27.8019单位后击中B船
% 导弹飞行的时间为2.7802分钟
% 击中点的坐标：
%    26.5523    6.5523


function dy=df5(t,y)
    v = 200; 
    dy=zeros(2,1);
    dy(1)=3*v*(20+sqrt(2)/2*v*t-y(1))/sqrt((20+sqrt(2)/2*v*t-y(1))^2+(sqrt(2)/2*v*t-y(2))^2);
    dy(2)=3*v*(sqrt(2)/2*v*t-y(2))/sqrt((20+sqrt(2)/2*v*t-y(1))^2+(sqrt(2)/2*v*t-y(2))^2);
end

微分方程模型

image-20240129141545307
image-20240129141735253
image-20240129141815500
image-20240129142046760
image-20240129142146748
image-20240129142209940
image-20240129142557126
image-20240129142619881

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%% Malthus模型（马尔萨斯模型）
clear;clc
x = dsolve('Dx=r*x','x(0)=x0','t')    % x = dsolve('Dx=r*x','x(t0)=x0','t')
% x = x0*exp(r*t)
% 怎么把上面这个式子中的x0和r替换成确定的值？
x0 = 100; 
r = 0.1;
subs(x)

% 初始人口为1000，画出50年且增长率分别为0.5%，1%，1.5% 一直到5%的人口变化曲线
x0 = 1000;
for i = 1:10
    r = 0.005*i;
    xx = subs(x);
    fplot(xx,[0,50])   % fplot函数可以绘制表达式的图形
    hold on
end


%% 阻滞增长模型（logistic模型）
clear;clc
x = dsolve('Dx=r*(1-x/xm)*x','x(t0)=x0','t')   % 化简后和书上的结果一样
% mupad  % 把计算出来的结果粘贴过去可以得到直观的表达式
% 高版本可以使用实时脚本
t0 = 0;
x0 = 1000;
xm = 10000;
r = 0.05;
xx = subs(x);    %  10000/(exp(log(9) - t/20) + 1)
fplot(xx,[0,200])  

% 如果不会用上面的fplot函数怎么办？
t = 0:200;
x = 10000 ./ (exp(log(9) - t/20) + 1);
plot(t,x,'-')

预测美国人口

image-20240129143246601
image-20240129143301788 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
clear;clc
year = 1790:10:2000;
population = [3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,179.3,204.0,226.5,251.4,281.4];
cftool  % 拟合工具箱
% (1) X data 选择 year
% (2) Y data 选择 population
% (3) 拟合方式选择：Custom Equation (自定义方程)
% (4) 修改下方的方框为：x = f(t) = xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))
% (5) 左边的result一栏最上面显示：Fit computation did not converge:即没有找到收敛解，右边的拟合图形也表明拟合结果不理想
% (6) 点击Fit Options，修改非线性最小二乘估计法拟合的初始值(StartPoint), r修改为0.02，xm修改为500
% (7) 此时左边的result一览得到了拟合结果：r = 0.02735, xm = 342.4
% (8) 依次点击拟合工具箱的菜单栏最左边的文件—Generate Code(导出代码到时候可以放在你的论文附录)，可以得到一个未命名的脚本文件
% (9) 在这个打开的脚本中按快捷键Ctrl+S，将这个文件保存到当前文件夹。
% (10) 在现在这个文件中调用这个函数得到参数的拟合值和预测的效果
[fitresult, gof] = createFit(year, population) 
t = 2001:2030;
xm = 342.4;   
r =  0.02735;
predictions = xm./(1+(xm./3.9-1).*exp(-r.*(t-1790)));  % 计算预测值（注意这里要写成点乘和点除）
figure(2)
plot(year,population,'o',t,predictions,'.')  % 绘制预测结果图
disp(predictions)  % 预测的数值




function [fitresult, gof] = createFit(year, population)
%CREATEFIT(YEAR,POPULATION)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input : year
%      Y Output: population
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  另请参阅 FIT, CFIT, SFIT.

%  由 MATLAB 于 02-Jan-2020 23:14:10 自动生成


%% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( year, population );

% Set up fittype and options.
ft = fittype( 'xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))', 'independent', 't', 'dependent', 'x' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.2 500];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, xData, yData );
legend( h, 'population vs. year', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
xlabel year
ylabel population
grid on

捕食者猎物模型

背景

image-20240129143537340
image-20240129143946435
image-20240129145756075
image-20240129145847371
image-20240129145919864
image-20240129145944577 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
clear;clc
%  Matlab求不出来解析解
% dsolve('Dx1=x1*(0.9-0.1*x2)','Dx2=x2*(-0.6+0.02*x1)','x1(0)=25,x2(0)=2','t')

% 下面用ode45函数求数值解
% 自变量t的范围为0-15年，食饵和捕食者（鲨鱼）初始值分别为25和2 
% 注意：战前和战后是战争开始前和战争开始后的简写哦
[t1,x1]=ode45('pre_war',[0 15],[25 2]);  % 战前的微分方程       
[t2,x2]=ode45('past_war',[0 15],[25 2]);     % 战后的微分方程
% [t1,x1]=ode45('pre_war',[0:15],[25 2]);  % 战前的微分方程       
% [t2,x2]=ode45('past_war',[0:15],[25 2]);     % 战后的微分方程

pre_prey=x1(:,1); pre_shark=x1(:,2);   % 战前的食饵和鲨鱼的数量
past_prey=x2(:,1); past_shark=x2(:,2);  % 战后的食饵和鲨鱼的数量

figure(1)
plot(pre_prey,pre_shark,'--r',past_prey,past_shark,'-b')
legend('战前','战后')
title('鲨鱼和食饵数量变化的相轨线图')
xlabel('食饵数量');   ylabel('鲨鱼数量')

figure(2)
plot(t1,pre_prey,'--r',t1,pre_shark,'-r',t2,past_prey,'--b',t2,past_shark,'-b')
legend('战前食饵数量','战前鲨鱼数量','战后食饵数量','战后鲨鱼数量')
xlabel('时间');   ylabel('数量')

% 鲨鱼比例  =  鲨鱼数量 /（鲨鱼数+食饵数）
pre_rate=pre_shark./(pre_prey+pre_shark);  % 战前的鲨鱼比例    
past_rate=past_shark./(past_prey+past_shark);   % 战后的鲨鱼比例
figure(3)
plot(t1,pre_rate,'--r',t2,past_rate,'-b')
legend('战前的鲨鱼比例','战后的鲨鱼比例')





% 战后
function dx=past_war(t,x)
    dx=zeros(2,1); 
    dx(1)=x(1)*(0.9-0.1*x(2));
    dx(2)=x(2)*(-0.6+0.02*x(1));
end


% սǰ
function dx=pre_war(t,x)
    dx=zeros(2,1); 
    dx(1)=x(1)*(0.7-0.1*x(2));
    dx(2)=x(2)*(-0.8+0.02*x(1));
end

种群相互竞争模型

image-20240129150258790
image-20240129150851048
image-20240129151411054
image-20240129151639245
image-20240129151742576 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
clc;clear
%  Matlab求不出来解析解
% dsolve('Dx1 = 0.5*x1*(1-x1/300-0.5*x2/500)','Dx2=0.5*x2*(1-x2/500-2*x1/300)','x1(0)=80,x2(0)=100','t')

% 下面用ode45函数求数值解
% 自变量为时间t，范围为0-30； 甲乙两个种群的数量初始值为80，100（随便给的，大家可以调整来看结果的变化）
[t,x]=ode45('fun',[0 30],[80 100]); 
plot(t,x(:,1),'r-',t,x(:,2),'b-')  % x的第一列是甲种群数量，x的第二列是乙种群数量
legend('种群甲','种群乙')
% axis([0 30 0 500])


function dx=fun(t,x)   % 大家可以修改里面的参数，来看结果的变化
    r1=0.5;  r 2=0.5; % 甲乙的增长率
%     r1=0.8;  r2=1; % 甲乙的增长率
    N1=300;   N2=500;   % 甲乙的最大数量
    % sigma1: 单位数量的乙种群(相对于N2)消耗的供养甲的食物量为单位数量的甲(相对于N1)消耗的供养甲的食物量的倍数。
    % sigma2: 单位数量的甲种群(相对于N1)消耗的供养乙的食物量为单位数量的乙(相对于N2)消耗的供养乙的食物量的倍数。
    sigma1=0.5;  sigma2=2;   
%     sigma1=0.5;  sigma2=4;   
%     sigma1=0.4;  sigma2=0.2;
%   当sigma1和sigma2同时大于1时（这种现象本身在自然界就几乎不可能出现），得到的结果不稳定。
%     sigma1=3;  sigma2=2;   
%     sigma1=2.2;  sigma2=2;   

    dx = zeros(2,1);
    dx(1) = r1*x(1)*(1-x(1)/N1-sigma1*x(2)/N2);
    dx(2) = r2*x(2)*(1-x(2)/N2-sigma2*x(1)/N1);
end

种群相互依存模型

image-20240129152337478
image-20240129152352501

情况一：

image-20240129152422993

情况二：

image-20240129152504494

情况三：

image-20240129152541262

情况一图像：

image-20240129152609368
image-20240129152614946
image-20240129152744579
image-20240129152809888

情况二图像：

image-20240129152835035

情况三图像：

image-20240129153016550

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
clc;clear
% 情况一：甲可以独自生存，乙不能独自生存
[t,x]=ode45('fun1',[0 50],[80 100]); 
figure(1)
plot(t,x(:,1),'r-',t,x(:,2),'b-')  % x的第一列是甲种群数量，x的第二列是乙种群数量
legend('种群甲','种群乙')

% 情况二：甲乙均可以独自生存
[t,x]=ode45('fun2',[0 50],[80 100]); 
figure(2)
plot(t,x(:,1),'r-',t,x(:,2),'b-')  % x的第一列是甲种群数量，x的第二列是乙种群数量
legend('种群甲','种群乙')

% 情况三：甲乙均不能独自生存
[t,x]=ode45('fun3',[0 50],[80 100]); 
figure(3)
plot(t,x(:,1),'r-',t,x(:,2),'b-')  % x的第一列是甲种群数量，x的第二列是乙种群数量
legend('种群甲','种群乙')






% 情况一：甲可以独自生存，乙不能独自生存
function dx=fun1(t,x)   % 大家可以修改里面的参数，来看结果的变化
    r1=0.5;  r2=0.5; % 甲的增长率和乙的死亡率
    N1=300;   N2=500;   % 甲乙的最大数量
    % sigma1: 单位数量的乙种群(相对于N2)提供的供养甲的食物量为单位数量的甲(相对于N1)消耗的供养甲的食物量的倍数。
    % sigma2: 单位数量的甲种群(相对于N1)提供的供养乙的食物量为单位数量的乙(相对于N2)消耗的供养乙的食物量的倍数。
    sigma1=0.2;  sigma2=2;   
%     sigma1=0.2;  sigma2=0.8;
% 注意：当sigma1*sigma2>1时，微分方程不稳定，matlab计算数值解时可能会报错，这时候需要调整计算的范围。
%     sigma1=3;  sigma2=3;   
    dx = zeros(2,1);
    dx(1) = r1*x(1)*(1-x(1)/N1+sigma1*x(2)/N2);
    dx(2) = r2*x(2)*(-1-x(2)/N2+sigma2*x(1)/N1);
end





% 情况二：甲乙均可以独自生存
function dx=fun2(t,x)   % 大家可以修改里面的参数，来看结果的变化
    r1=0.5;  r2=0.5; % 甲的增长率和乙的增长率
    N1=300;   N2=500;   % 甲乙的最大数量
    % sigma1: 单位数量的乙种群(相对于N2)提供的供养甲的食物量为单位数量的甲(相对于N1)消耗的供养甲的食物量的倍数。
    % sigma2: 单位数量的甲种群(相对于N1)提供的供养乙的食物量为单位数量的乙(相对于N2)消耗的供养乙的食物量的倍数。
    sigma1=0.2;  sigma2=2;   
%     sigma1=0.2;  sigma2=0.8;
    % 注意：当sigma1*sigma2>1时，微分方程不稳定，matlab计算数值解时可能会报错。
    dx = zeros(2,1);
    dx(1) = r1*x(1)*(1-x(1)/N1+sigma1*x(2)/N2);
    dx(2) = r2*x(2)*(1-x(2)/N2+sigma2*x(1)/N1);
end


% 情况三：甲乙均不能独自生存
function dx=fun3(t,x)   % 大家可以修改里面的参数，来看结果的变化
    r1=0.2;  r2=0.2; % 甲的死亡率和乙的死亡率
    N1=300;   N2=500;    % 甲乙的最大数量
    % sigma1: 单位数量的乙种群(相对于N2)提供的供养甲的食物量为单位数量的甲(相对于N1)消耗的供养甲的食物量的倍数。
    % sigma2: 单位数量的甲种群(相对于N1)提供的供养乙的食物量为单位数量的乙(相对于N2)消耗的供养乙的食物量的倍数。
    sigma1=0.2;  sigma2=2;   
%     sigma1=5;  sigma2=5; 
%     sigma1=10;  sigma2=10;  % 这时候甲乙两个种群都能活下去了 
    dx = zeros(2,1);
    dx(1) = r1*x(1)*(-1-x(1)/N1+sigma1*x(2)/N2);
    dx(2) = r2*x(2)*(-1-x(2)/N2+sigma2*x(1)/N1);
end