目录

Matlab求解微分方程

解析解

数值解

一阶微分方程

高阶微分方程

人口预测模型

Malthus 马尔萨斯模型

Logistic 阻滞增长模型

matlab曲线拟合工具箱 Curve fitting

捕食者猎物模型

结果分析

种群竞争模型

种群依存模型

传染病模型

SI模型

参数β(接触率)降低的情况

考虑出生率与死亡率

只考虑疾病的死亡率

​编辑

考虑出生率、死亡率、因病死亡率

SIS模型

康复率提升

出生率、死亡率、因病死亡率

SIR模型

康复率增加

因病死亡率

SIRS模型

因病死亡

完全和不完全康复

SEIR模型

潜伏者没有传染性

潜伏者有传染性


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 解

结果分析

相轨线图:

(1) 从下方开始逆时针去看上面这个相轨线图 :随着食饵增多,则鲨鱼易于取食,所以鲨鱼数量增加;然而,由于鲨鱼数量的增多而需要食用更多的食饵,此时食饵的数量正在下降,直到饵数量特别少时,鲨鱼进入了饥饿的状态而使得鲨鱼的数量急剧下降,这时部分食饵得以存活,食饵数量重新开始回升;食饵与鲨鱼的数量交替增减形成了生物圈中的动态平衡。
(2) 另外,蓝色的圈代表战后的数量变化,可以看出,蓝色的圈几乎把红色的圈包在里面,且蓝色圈的上方和红色圈的上方差距较大,这说明战后鲨鱼的增长速度要明显快于红圈,这导致了鲨鱼数量的大幅增长,和生物学家 D'Ancona 观察到的现象相一致。

种群竞争模型

竞争存在于资源有限的情况下,所以基于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 > 60beta = beta/10; % 第60期后禁止大规模聚会,使得传染强度beta缩小为原来的10倍end
%     beta = 0.1 - 0.001*t
%     if beta < 0
%         beta = 0;
%     enddx = zeros(2,1); % x(1)表示S  x(2)表示Idx(1) = - beta*x(1)*x(2)/TOTAL_N;  dx(2) = beta*x(1)*x(2)/TOTAL_N;
end

考虑出生率与死亡率

只考虑疾病的死亡率

考虑出生率、死亡率、因病死亡率

SIS模型

治愈后又再患病,治愈率为 α

康复率提升

医疗设备升级、建立医院

出生率、死亡率、因病死亡率

与SI类似

SIR模型

这里的康复指的是产生 抗体,不会再感染,前面的SIS是没有产生抗体的暂时恢复

康复率增加

疫苗、医疗设备升级

因病死亡率

SIRS模型

因病死亡

完全和不完全康复

SEIR模型

潜伏者没有传染性

潜伏者有传染性

把未感染者被感染的概率分为 被潜伏者感染、被感染者感染

数学建模:17 微分方程相关推荐

  1. 求微分方程的平衡点matlab,数学建模之微分方程建模与平衡点理论

    微分方程 列微分方程常用的方法: (1)根据规律列方程 利用数学.力学.物理.化学等学科中的定理或经过实验检验的规律来建立微分方程模型. (2)微元分析法 利用已知的定理与规律寻找微元之间的关系式,与 ...

  2. 数学建模之微分方程(符实现例题和MATLAB源码)

    微分方程的基本概念 微分方程:一般的,凡表示未知函数.未知函数的导数与自变量之间的关系的方程,叫做微分方程,有时也简称方程. 微分方程的阶:微分方程中所出现的未知函数的最高阶导数的阶数,叫做微分方程的 ...

  3. 数学建模微分方程导弹问题matlab求解,数学建模之微分方程(符实现例题和MATLAB源码)...

    微分方程的基本概念 微分方程:一般的,凡表示未知函数.未知函数的导数与自变量之间的关系的方程,叫做微分方程,有时也简称方程. 微分方程的阶:微分方程中所出现的未知函数的最高阶导数的阶数,叫做微分方程的 ...

  4. 数学建模:微分方程模型—常微分方程数值解算法及 Python 实现

    目录 一.显式欧拉 (Euler) 法 二.显式欧拉法的改进 隐式欧拉法 (后退欧拉法) 梯形法 两步欧拉法 (中点法) 预报 - 校正法 (改进欧拉法) 三.龙格 - 库塔 (Runge-Kutta ...

  5. 数学建模【微分方程模型(介绍、分析方法、数值模拟、传染病问题的建模和分析、经济增长模型、人口增长预测和控制模型)】

  6. 2019年上海市数学建模讲座(3)微分方程建模方法

    第三场微分方程建模讲座笔记 主讲人:董程栋,上海财经大学数学学院 微分方程: 定义:联系着自变量,未知函数与它的导数之间的关系式 物体冷却过程中的数学模型 牛顿冷却定律:物体温度变化速度与物体和介质温 ...

  7. Python小白的数学建模课-17.条件最短路径算法

    条件最短路径问题,指带有约束条件.限制条件的最短路径问题.例如: 顶点约束,包括必经点或禁止点的限制: 边的约束,包括必经路段.禁行路段和单向路段:无权路径长度的限制,如要求经过几步或不超过几步到达终 ...

  8. Python小白的数学建模课-10.微分方程边值问题

    小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型边值问题的建模与求解,不涉及算法推导和编程,只探讨如何使用 Pytho ...

  9. Python小白的数学建模课-09.微分方程模型

    小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型的建模与求解,通过常微分方程.常微分方程组.高阶常微分方程 3个案例手 ...

  10. 种群内禀增长率matlab求法,数学建模讲义:第三讲微分方程模型

    <数学建模讲义:第三讲微分方程模型>由会员分享,可在线阅读,更多相关<数学建模讲义:第三讲微分方程模型(74页珍藏版)>请在人人文库网上搜索. 1.第三讲 微分方程模型,动态模 ...

最新文章

  1. 为啥我从后台查到的值在页面显示的是undefined_再谈一个管理后台列表功能应有的素质...
  2. Java调用PHP,跑PHP代码
  3. SAP MM PIR里的Lower Limit Upper Limit
  4. python上海培训哪里比较好-上海python培训哪家好 Python需要多久学会
  5. 极简数据分析实操指南(上)
  6. 【数值分析】顺序高斯消去法和列主元高斯消去法的三个主要不同点
  7. MFC线程自定义消息
  8. 3.1 Ext JS 组件总览
  9. React.js 小书 Lesson1-2 - 前端组件化(一):从一个简单的例子讲起
  10. invalid use of constructor as a template 编译错误
  11. VoLTE业务端到端流程
  12. IOS开发之UI进阶(设置圆角,边框颜色,边框宽度)
  13. 转:make cmake和catkin_make的区别
  14. segment:?co?_如何跟踪用户动作并了解它们:Segment + MixPanel
  15. git diff:Linux使用meld做git的diff工具
  16. 夜神模拟器连接Android Studio
  17. QGraphicsItem图元拖动绘制(二)
  18. J2EE三层架构简介
  19. 【大学生数学竞赛】公式大全(补充中)
  20. python之红楼梦词频统计并生成图云

热门文章

  1. 拉别人或者公司的Uni-app项目,在HBuilder运行小程序时打开了微信开发者工具,但微信开发者工具未运行项目
  2. 无线显示android,手机变曲面无线显示
  3. 【C++】字符串详解
  4. 获取iPhone本机IP地址新方法
  5. unix、Linux知多少
  6. Spotify Music Converter for Mac如何注册?
  7. 双击bat文件变成打开文本编辑器该如何解决
  8. SAP 各种错误解决大全
  9. void xiaob::printwhat(string* str) 如果括号是* 外面调用传递什么,如果括号是 外面调用传递什么...
  10. 飞桨开源深度学习平台介绍