该博客为个人学习清风建模的学习笔记,代码全部摘自清风老师,部分课程可以在B站:【强烈推荐】清风:数学建模算法、编程和写作培训的视频课程以及Matlab等软件教学_哔哩哔哩_bilibili
目录
名称 | 重要性 | 难度 |
微分方程模型 | ★★★★★ | ★★★★ |
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模型