清风数学建模——模型学习微分方程模型

微分方程及微分方程模型的基础概念

Matlab求微分方程的解析解

Matlab求一阶微分方程的数值解

Matlab求高阶微分方程的数值解

例题2：刚性问题

例题3：导弹追击问题

微分方程模型之：人口预测模型

微分方程模型之：捕食者-猎物模型

微分方程模型之：人口预测模型微分方程模型之：种群相互竞争模型、相互依存模型

微分方程模型之：传染病SI模型

非原创，看清风数学建模视频课做的笔记，建议去B站观看清风数学建模的课程或者购买正课~

基础概念：
常微分方程：y’+4y-2=0；导数的方程，一个变量
偏微分方程：偏导数的方程，多个变量

Matlab求微分方程的解析解：dsolve(‘方程1’,‘方程2’，…，‘方程n’,‘初始条件’,‘自变量’)；一阶微分方程；含未知系数的微分方程；高阶微分方程及初始条件；微分方程组

1、Matlab在表示微分方程时，用字母D表示微分，D2,D3分别表示二阶、三阶微分，后面跟的是要求解的因变量；
2、必须要指定自变量，要不然Matlab默认自变量是t，若实际自变量不是t，则计算结果是不正确的；
3、初始值可以不给，不给求出来的就是通解
4、乘号* 不能省略
%% 求解y-y’=2x
clear;clc
dsolve(‘y-Dy=2*x’,‘x’)
% 新版本的Matlab可能会淘汰上面的这种写法，但目前就用这种写法吧

%% 求解含有未知系数的微分方程 y-y’=ax
clear;clc
dsolve(‘y-Dy=a*x’,‘x’)

%% 求解二阶微分方程，已知初始条件 y’‘+4y’+29y=0 且 y(0)=0 y’(0)=15
clear;clc
dsolve(‘D2y+4Dy+29y=0’,‘y(0)=0,Dy(0)=15’,‘x’)

%% 求解微分方程组

clear;clc
[x,y,z]=dsolve(‘Dx=2x-3y+3z+t’,'Dy=4x-5y+3z+t’,‘Dz=4x-4y+2*z+t’,‘t’)

Matlab求一阶微分方程的数值解：
微分方程的数值解函数如下： [x,y]=solver(‘f’,ts,x0,options)
1、x代表自变量； y代表函数值； solver代表求解函数，总共有七种，如下图所示。ts=[t0,tf]代表自变量的初值和终值；x0表示函数的初始值； f是待求解的微分方程需要用函数文件编写，且要写成标准形式，如下图所示
并用@进行调用，或者用单引号；options = odeset(‘reltol’,1e-4,‘abstol’,1e-8) 用来设置相对误差和绝对误差，相对误差是reltol，绝对误差是abstol
2、使用Matlab求解数值解时，高阶微分方程组必须转换为一阶微分方程组；
3、优先使用ode45;解不出来用ode15s.

刚性：短时间内数值变化很大，如下所示

非刚性：如下所示

实例如下：

%% 微分方程的数值解函数如下： [x,y]=solver(‘f’,ts,x0,options)

%% 求方程y-y’=2x的数值解，初始值y(0)=3, x范围是[0,2]
%标准形式为：y’=y-2x; 将其写为df1.m文件；
% function dy=fun1(x,y)
% dy=y-2x;
% end
clear;clc
[x,y]=ode45(‘fun1’,[0,2],3)
figure(1)
plot(x,y,'r’)
% 与解析解进行对比
dsolve(‘y-Dy=2*x’,‘y(0)=3’,‘x’)
xx=0:0.1:2
yy=2.*xx + exp(xx) + 2;
hold on
plot(xx,yy,‘b-’)%解析解

%% 求解微分方程组的数值解

% function dy=fun2(x,y)
% dy=zeros(3,1);
% dy(1)=y(2)*y(3);
% dy(2)=-y(1)y(3);
% dy(3)=-0.51y(1)y(2);
% end
[x,y]=ode45(‘fun2’,[0,4pi],[0,1,1])
%y有三列，第一列是y1,第2列是y2,第3列是y3

Matlab求高阶微分方程的数值解：

例题1如下：

需要将高阶微分方程转为一阶微分方程，降阶的方法为变量替换，即令y1=y,y2=y’，如下所示，降阶为一阶微分方程组，降阶的过程只能自己手算；

%% 求高阶微分方程的数值解
% 例题1
% 先看有没有解析解
clear;clc
dsolve(‘(1+x^2)D2y=2xDy’,‘y(-2)=3,Dy(-2)=4’,‘x’)
% y=(4x*(x^2 + 3))/15 + 101/15;
x=-2:0.1:2;
y=(4.x.(x.^2 + 3))./15 + 101/15;
plot(x,y,‘r’)
% 再看有没有数值解
[xx,yy]=ode45(‘fun3’,[-2,2],[3,4])
hold on
plot(xx,yy(:,1),‘b*’)% yy的第一列为待求解的y,第二列为y’,不关心；

例题3：导弹追击问题

1、相撞时，导弹与船的距离d=0，但对于数值法求解时，应设置d≈0’可以设置为小于0.001或者更小
%导弹轨迹
clear;clc
options = odeset(‘reltol’,1e-4,‘abstol’,1e-8);
[t,y]=ode45(‘fun_daodan’,[0,0.1],[0,0],options);
plot(y(:,1),y(:,2))
%船的轨迹
hold on
plot([20,30],[0,10])
%判断相撞
x=y(:,1);
y=y(:,2);
n=length(t);
d=0;
s=0;
for i=1:n
d(i)=abs(x(i)-y(i)-20)./sqrt(2);
if d(i)<0.001
for j=1:i-1
s=s+sqrt((x(j+1)-x(j))^{2+(y(j+1)-y(j))}2);
end
disp([‘导弹运行’,num2str(s),‘单位后与船撞击’])
disp([‘导弹运行时间为’,num2str(t(i)60)])
break
end
end
hold on
plot(x(i),y(i),'')
text(x(i),y(i),‘击中点’)

微分方程模型之：人口预测模型
%% 人口预测模型
% Malthus模型
%求解含有未知量的微分方程的解析解，不需要设置为符号变量，直接写方程即可
clear;clc
dsolve(‘Dx=rx’,‘x(0)=x0’,‘t’)
%求解得到：x(t)=x0exp(rt)
% Logistic模型
clear;clc
x=dsolve('Dx=r(1-x/xm)*x’,‘x(t0)=x0’,‘t’)
simplify(x)%化简表达式
t0 = 0;
x0 = 1000;
xm = 10000;
r = 0.05;
xx=subs(x)%将未知量替换为具体的数值
fplot(xx,[0,200])

%% 例题：美国人口预测

clear;clc
% 思路：先求出微分方程模型的解析式，然后使用cftool的自定义方程拟合得到解析式中的未知量，最后进行预测
t=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];
x=dsolve(‘Dx=r*(1-x/xm)*x’,‘x(1790)=3.9’,‘t’);
xx=simplify(x);
% cftool
[fitresult, gof] = createFit(t, population)
r=0.02735;
xm=342.4;
tp=2001:1:2030;
pre=(39.*xm.exp(r.(tp - 1790)))./(10.*xm + 39.exp(r.(tp - 1790)) - 39)
plot(t,population,‘r’,tp,pre,‘b’)

微分方程模型之：捕食者-猎物模型
本质还是求微分方程，然后画图
%% 捕食者-猎物模型
% 求解战前和战后鲨鱼与猎物的数量
clear;clc
[t1,x1]=ode45(‘pre_war’,[0 15],[25,2]);
[t2,x2]=ode45(‘post_war’,[0 15],[25,2]);
pre_liewu=x1(:,1);
pre_shark=x1(:,2);
post_liewu=x2(:,1);
post_shark=x2(:,2);
% 相轨图
figure(1)
plot(pre_liewu,pre_shark,‘–r’,post_liewu,post_shark,‘-b’)
legend(‘战前’,‘战后’)
% 鲨鱼猎物数量
figure(2)
plot(t1,pre_liewu,‘–r’,t1,pre_shark,‘-r’)%战前
hold on
plot(t2,post_liewu,‘–b’,t2,post_shark,‘-b’);
legend(‘战前猎物数量’,‘战前鲨鱼数量’,‘战后猎物数量’,‘战后鲨鱼数量’)
% 鲨鱼比例
ratio_pre_shark=pre_shark./(pre_shark+pre_liewu);
ratio_post_shark=post_shark./(post_shark+post_liewu);
figure(3)
plot(t1,ratio_pre_shark,‘–r’)
hold on
plot(t2,ratio_post_shark,‘-b’)
legend(‘战前鲨鱼比例’,‘战后鲨鱼比例’)

微分方程模型之：种群相互竞争模型、相互依存模型
%% 种群相互竞争模型
% 看下有没有解析解
% clear;clc
% r1=0.5;
% r2=0.5;
% N1=300;
% N2=500;
% sigma1=0.5;
% sigma2=2;
% [x1,x2]=dsolve(‘Dx1=r1x1(1-x1/N1-sigma1x2/N2)', …
% 'Dx2=r2x2*(1-x2/N2-sigma2*x1/N1)’, …
% ‘x1(0)=80,x2(0)=100’,‘t’)
% 求不出来解析解
% 数值解
clear;clc
[t,x]=ode45(‘fun5’,[0,30],[80,100])
figure(1)
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’)
legend(‘种群甲’,‘种群乙’)

%% 种群相互依存模型
% 情况一：甲可以独立生存，乙不能独立生存
function dx=qingkuang1(t,x)
dx=zeros(2,1);
r1=0.5;
r2=0.5;
N1=300;
N2=500;
sigma1=0.2;
sigma2=2;
dx(1)=r1.x(1).(1-x(1)./N1+sigma1.*x(2)./N2);
dx(2)=r2.x(2).(-1+sigma2.*x(1)./N1-x(2)./N2);
end
clear;clc
[t,x]=ode45(‘qingkuang1’,[0,50],[80,100])
figure(1)
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’)
legend(‘种群甲’,‘种群乙’)
% 情况二：甲可以独立生存，乙可以独立生存
clear;clc
[t,x]=ode45(‘qingkuang2’,[0,50],[80,100])
figure(1)
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’)
legend(‘种群甲’,‘种群乙’)
% 情况三：均不能独立生存
clear;clc
[t,x]=ode45(‘qingkuang3’,[0,50],[80,100])
figure(1)
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’)
legend(‘种群甲’,‘种群乙’)

微分方程模型之：传染病SI模型
%% 传染病SI模型
%% 改进1：考虑第60期后禁止大规模聚会，传染强度beta缩小为原来的10倍
clear;clc
function dx=SI_model(t,x)
dx=zeros(2,1)
global N beta
dx(1)=-beta.*x(1).*x(2)./N;
dx(2)=beta.*x(1).*x(2)./N;
end

global N beta
N=1000;% 总人数
beta=0.1;% 初期传染强度
i0=1; %初始感染者人数
[t,x]=ode45(‘SI_model’,[0,60],[N-i0,i0])% 初期计算数据
beta=beta/10;%后期传染强度
post_S=x(end,1); % 后期初始S
post_I=x(end,2); % 后期初始I
[post_t,post_x]=ode45(‘SI_model’,[61,200],[post_S,post_I]) % 后期计算数据
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’)
hold on
plot(post_t,post_x(:,1),‘-r’,post_t,post_x(:,2),‘-b’)
% 上述方法为方法1，方法2为，再微分方程SI_model中加入if判断句，若if>60，则beta=beta/10;
% 而后只需要[t,x]=ode45(‘SI_model’,[0,20],[N-i0,i0]) 即可

%% 改进2：增加人口自然出生率和死亡率，但不考虑疾病的死亡率
function dx=SI_model2(t,x)
dx=zeros(3,1);
global beta N u v
% N=x(1)+x(2);
dx(1)=-beta.*x(1).x(2)./(x(1)+x(2))+u.(x(1)+x(2))-v.*x(1);
dx(2)=beta.*x(1).*x(2)./(x(1)+x(2))-v.*x(2);
dx(3)=v.*x(1)+v.*x(2);
end
% 也可以写成下面这样
% function dx=SI_model2(t,x)
% dx=zeros(3,1);
% global beta N u v
% N=x(1)+x(2);
% dx(1)=-beta.*x(1).*x(2)./N+u.*N-v.*x(1);
% dx(2)=beta.*x(1).*x(2)./N-v.*x(2);
% dx(3)=v.*x(1)+v.*x(2);
% end
clear;clc
global beta N u v
beta=0.1;
u=0.002;
v=0.001;
N=1000;
i0=1;
[t,x]=ode45(‘SI_model2’,[0,200],[N-i0,i0,0])
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’,t,x(:,3),‘-k’)

%% 改进3：不考虑人口自然出生率和死亡率，只考虑疾病的死亡率
function dx=SI_model3(t,x)
dx=zeros(3,1);
global beta N d
N=x(1)+x(2);
dx(1)=-beta.*x(1).*x(2)./N;
dx(2)=beta.*x(1).*x(2)./N-d.*x(2);
dx(3)=d.*x(2);
end
clear;clc
global beta N d
N=1000;
i0=1;
beta=0.1;
d=0.01;
[t,x]=ode45(‘SI_model3’,[0,500],[N-i0,i0,0])
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’,t,x(:,3),‘-k’)

%% 改进4：同时考虑人口自然出生率和死亡率和疾病的死亡率.
function dx=SI_model4(t,x)
dx=zeros(4,1);
global beta N d u v
N=x(1)+x(2);
dx(1)=-beta.*x(1).*x(2)./N+u.*N-v.*x(1);
dx(2)=beta.*x(1).*x(2)./N-d.*x(2)-v.*x(2);
dx(3)=d.*x(2);
dx(4)=v.*x(1)+v.*x(2);
end
clear;clc
global beta N d u v
N=1000;
i0=1;
d=0.01;
beta=0.1;
u=0.002;
v=0.001;
[t,x]=ode45(‘SI_model4’,[0,500],[N-i0,i0,0,0])
plot(t,x(:,1),‘-r’,t,x(:,2),‘-b’,t,x(:,3),‘-k’,t,x(:,4),‘-g’)