微分方程模型★★★★★

该博客为个人学习清风建模的学习笔记,代码全部摘自清风老师,部分课程可以在B站:【强烈推荐】清风:数学建模算法、编程和写作培训的视频课程以及Matlab等软件教学_哔哩哔哩_bilibili

目录

1引例

2概述

2.1基本概念

2.2如何建立微分方程

3MATLAB求解解析解

3.1dsolve函数

3.2例题

4MATLAB求解数值解

4.1solver函数

4.2例题

5人口预测模型

5.1马尔萨斯模型

5.2阻滞增长模型

5.3例题

6捕食者-猎物模型

7种群相互竞争模型

8种群相互依次模型

9传染病模型

9.1SI模型

9.2SIS模型

9.3SIR模型

9.4SRIS模型

9.5SEIR模型


名称重要性难度
微分方程模型★★★★★★★★★

1引例

 位于坐标原点的A船向位于其正东⽅20个单位的B 船发射导弹,导弹始终对准B船,B船以时

V 单位(常数)沿东北⽅向逃逸。若导弹的速度为 3V ,导弹的射程是 50 个单位,画出导弹运⾏的曲线,导弹是否能在射程内击中B 船?

2概述

2.1基本概念

2.2如何建立微分方程

建立模型:专业知识+套用现有的模型

3MATLAB求解解析解

3.1dsolve函数

dsolve('a,'b','x')

a:方程

b:初始条件

x:自变量

3.2例题

%% 例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)

%% 例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)

4MATLAB求解数值解

4.1solver函数

[x,y] = solver('f',[x0,xn],y0,option)

x:自变量

y:函数值

solver:求解函数

f:待解的微分方程

x0:x的初值

xn:x的终值

y0:y的初值

option:精度设置(相对误差和绝对误差) 

4.2例题

%% 调用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);
figure(1)
plot(x,y,'b-') 


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

%% 在真正比赛中,我们往往会倾向于先去看有无解析解,如果没有解析解再考虑数值解
% [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 :0.0001: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

其他例题的代码就不放在博客了,逻辑写法差不多,会使用solver函数就可以解出来。注意在解微分方程的时候先看一下是否可以求出解析解,如果不行再解数值解。 

下面是引例的代码,摘自清风老师:

clear; clc
v = 200;
[t,y]=ode45('df5' ,[0,0.1],[0 0]); 
x = y(:,1);  y =y(:,2);
plot(x, y,'*','MarkerSize',1)  
hold on
plot([20,30],[0,10]) 
hold on
% 下面我们找到导弹与B船相撞的点
n =length(t);  
for i = 1:n 
    d = t(i) * v * 3;  % 计算这个时候导弹的飞行的距离
    % 下面我们直接计算导弹和B船的距离dd
    % 在t(i)这个时间点,导弹坐标为(x(i), y(i))B船坐标为(20+sqrt(2)/2*v*t(i), sqrt(2)/2*v*t(i))
    dd =sqrt((x(i)-(20+sqrt(2)/2*v*t(i)))^2+(y(i)-sqrt(2)/2*v*t(i))^2);  
    if dd < 0.01 % 如果这个距离足够小了,我们就认为相撞了
        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.01 % 如果射程大于50或导弹与B船的距离始终都没有小于0.01(这个数需要根据实际情况调整)
    disp('导弹没有击中B船');
end


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

5人口预测模型

5.1马尔萨斯模型

%% 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

5.2阻滞增长模型

%% 阻滞增长模型(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,'-')

5.3例题

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)  % 预测的数值

6捕食者-猎物模型

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]);     % 战后的微分方程

% 战前
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

% 战后
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

 7种群相互竞争模型

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('种群甲','种群乙')

function dx=fun(t,x)   % 大家可以修改里面的参数,来看结果的变化
    r1=0.5;  r2=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

8种群相互依次模型

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

9传染病模型

 9.1SI模型

 

9.2SIS模型

 

9.3SIR模型

9.4SRIS模型

 

9.5SEIR模型

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值