目录
Matlab求解微分方程
先尝试求解析解(即函数表达式),求不出再求数值解(求出很多个数值,就能画图由离散的值近似原本的函数表达式图像)
解析解
% 方法1
dsolve('y-Dy=2*x','y(0)=3','x')
% 解析解:y = 2*x + exp(x) + 2
% Dy表示一阶微分 dy/dx
% D2y表示二阶微分
% 方法2
syms y(x)
eqn = (y - diff(y,x) == 2*x); % 要用双等号!!
cond = (y(0) == 3); % 初始条件
dsolve(eqn,cond)
% diff(y,x) 即 dy/dx
% 方法1
dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
% 3*sin(5*x)*exp(-2*x)
%%%%%%%% 画图:点乘!!!
x = -2:0.01:2;
y = 3*sin(5.*x)*exp(-2*x)
plot(x,y,'b-')
% 方法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)
% 方法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')
A = 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')
A.x
% 方法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)
求出解析解后,画图:记得把解析解中的加减乘除 换成 点乘、点乘方!!!
simplify(y) % simplify函数可以简化表达式
latex(y) % 转换成latex代码,复制到Axmath或者word自带的公式编辑器
数值解
[x,y] = ode45('df1',[0,2],3);
% [x,y] = ode45(@df1,[0,2],3); % 解决非刚性问题,这个解不出就试一试ode15s
% [x,y] = ode23('df1',[0,2],3);
% [x,y] = ode113('df1',[0,2],3);
% [x,y] = ode15s('df1',[0,2],3); % 解决刚性问题
%% 设定相对误差reltol和绝对误差,这样可以提高微分方程数值解的精度
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);
%%%%% 有参数的话只能用函数句柄传参
[t,x]=ode45(@(t,x) newfun2(t,x,GAMMA),[1:500],[N-i0 i0 0]);
一阶微分方程
先看有没有解析解,没有才求数值解
注意 .m 文件中要将微分方程写成标准形式(等式左边是一阶微分,右边是剩余部分)
%% 先看有无解析解,如果没有解析解再考虑数值解
[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.m,
% 自变量x的范围为[0,4*pi] ;初始值: y1(0)=0,y2(0)=y3(0)=1
[x, y] = ode45('df2', [0 4*pi], [0 1 1]);
% 返回值y是一个有3列的矩阵,第一列是y1、以此类推
% x 是范围中的一些取值,y是x对应的函数值,画出图像后用离散的数值解近似解析解的图像
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
%%% df2.m文件:
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
高阶微分方程
高阶微分方程求数值解,必须先转化为一阶微分方程(高数中学过的方法降阶)
% 先看这个微分方程有没有解析解
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的一阶导数)
%%%%%%%%%%%%%% df4.m
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
人口预测模型
Malthus 马尔萨斯模型
模型假设与推导:(绿框中为结果)
存在的问题:r 不是常数,随人口增大而减小,先假设是线性的递减函数(工程师原则)—— logistic模型
Logistic 阻滞增长模型
f(t) = xm/(1+(xm/3.9-1)*exp(-r*(t-1790)))
matlab曲线拟合工具箱 Curve fitting
自定义拟合方程:
在高级选择项中修改参数使得更加拟合:(修改了右下方的r、xm的起点、上下限)
捕食者猎物模型
只考虑一对捕食者与食饵:(基于指数增长的马尔萨斯模型)
食饵单独存在时候的增长率(受到捕食者数目的影响)、捕食者单独存在的死亡率(受到食饵数目的影响)、人工捕获食饵的影响(减少食饵增长率、增大捕食者死亡率)
食饵考虑增长率 dx/dt = rx 是正号,鲨鱼考虑死亡率 dx/dt = -rx 是负号!!!
食饵考虑被捕食的影响(负向影响),该影响与捕食者数目成正比,即 -λx;
捕食者考虑捕食的影响(正向影响),该影响与食饵数目成正比,即 +λx;
参数 x1 和 x2 的初值、r、e、λ 的值都自己确定,用 ode45 解
结果分析
相轨线图:
种群竞争模型
竞争存在于资源有限的情况下,所以基于logistic模型考虑种群自然生长的情况,还应考虑竞争者对自己食物的竞争能力(与对方的数目 x/N 成比例)
考虑的都是增长率
这里相乘 > 1 的话,画出来的图像不稳定!!
种群依存模型
考虑自身数目的阻滞作用(基于logistic模型)
传染病模型
易感者 S (Susceptible)
潜伏者 E (Exposed)
感染者 I (Infected)
康复者 R (Recovered)
基础模型:
总人数 N 只包括在传染病系统中的总人数,在SIR模型康复了不会再感染的人不包括在N中,所以在SIR模型中 N'=S+I 人数变化;
但在SI、SIS模型中 N = S+I 恒定;
但在SISRS模型康复了的人还会再患病,所以 N = S+I+R 恒定
SEIR中,由于有人完全康复所以人数变化 N'=S+E+I
扩展模型:
考虑出生率、死亡率、因病死亡率时,总人数不是恒定的,所以在SIR模型、SI、SIS模型中仍有 N'=S+I 人数变化;
SIRS中,考虑出生死亡因病死亡 N'=S+I;
考虑 R 分为完全康复 R1(不会再患病)、可能再患病 R2, N'=S+I+R2
规律:
- 有几个状态就有几个针对该状态的微分方程,方框中的是状态:对方框中的状态求微分 dS、dI、dE、dR、dND自然死亡、dID因病死亡
- 状态转移图中,输入线对应方程中的 + 号、输出线是 - 号
SI模型
不考虑人口的迁徙和出生死亡,人口总人数 N = S + I
微分方程:人为给定 β、N,求 S、I 的数值解并画出图像
参数β(接触率)降低的情况
function dx=fun2(t,x)
global TOTAL_N % 定义总人数为全局变量
beta = 0.1; % 易感染者与已感染者接触且被传染的强度
if t > 60
beta = beta/10; % 第60期后禁止大规模聚会,使得传染强度beta缩小为原来的10倍
end
% beta = 0.1 - 0.001*t
% if beta < 0
% beta = 0;
% end
dx = zeros(2,1);
% x(1)表示S x(2)表示I
dx(1) = - beta*x(1)*x(2)/TOTAL_N;
dx(2) = beta*x(1)*x(2)/TOTAL_N;
end
考虑出生率与死亡率
只考虑疾病的死亡率
考虑出生率、死亡率、因病死亡率
SIS模型
治愈后又再患病,治愈率为 α
康复率提升
医疗设备升级、建立医院
出生率、死亡率、因病死亡率
与SI类似
SIR模型
这里的康复指的是产生 抗体,不会再感染,前面的SIS是没有产生抗体的暂时恢复
康复率增加
疫苗、医疗设备升级
因病死亡率
SIRS模型
因病死亡
完全和不完全康复
SEIR模型
潜伏者没有传染性
潜伏者有传染性
把未感染者被感染的概率分为 被潜伏者感染、被感染者感染