数学建模-常微分方程模型及差分模型市公开课一等奖省名师优质课赛课一等奖课件

微分方程模型洛阳理工学院数理部1微分方程模型人口增加预测传染病模型种群模型2动态模型

描述对象特征随时间(空间)演变过程.

分析对象特征改变规律.

预报对象特征未来性态.

研究控制对象特征伎俩.

依据函数及其改变率之间关系确定函数.微分方程建模根据建模目和问题分析作出简化假设.

按照内在规律或用类比法建立微分方程.3对微分方程研究方法解在很广泛条件下存在，但能用有限解析式表示者极少.另辟它径：1、求数值解（近似解）；2、定性方法分析.4背景

年1625183019301960197419871999人口(亿)5102030405060世界人口增加概况中国人口增加概况

年1901192919531965198219901995人口(亿)4.265.486.027.2510.3211.3012.0012.95研究人口改变规律控制人口过快增加人口增加预测5指数增加模型——马尔萨斯提出(1798)惯用计算公式x(t)~时刻t人口基本假设

:人口(相对)增加率r

是常数今年人口x0,年增加率rk年后人口伴随时间增加，人口按指数规律无限增加6年17901800181018201830184018501860人口3.95.37.29.612.917.123.231.4年18701880189019001910192019301940人口38.650.262.976.092.0106.5123.2131.7年19501960197019801990人口150.7179.3204.0226.5251.4281.4-3-177Anna美国人口统计数据指数增加模型应用及不足

与19世纪以前欧洲一些地域人口统计数据吻合

适合用于19世纪后迁往加拿大欧洲移民后代

可用于短期人口增加预测

不符合19世纪后多数地域人口增加规律

不能预测较长久人口增加过程19世纪后人口数据人口增加率r不是常数(逐步下降)-3-178Anna阻滞增加模型(Logistic模型)人口增加到一定数量后，增加率下降原因：资源、环境等原因对人口增加阻滞作用且阻滞作用随人口数量增加而变大假设r~固有增加率(x很小时)xm~人口容量（资源、环境能容纳最大数量）r是x减函数-3-179Annadx/dtx0xmxm/2xmtx0x(t)~S形曲线,x增加先快后慢x0xm/2阻滞增加模型(Logistic模型)-3-1710Anna参数预计用指数增加模型或阻滞增加模型作人口预报，必须先预计模型参数r或r,xm

利用统计数据用最小二乘法作拟合例：美国人口数据（单位~百万）186018701880……196019701980199031.438.650.2……179.3204.0226.5251.4阻滞增加模型(Logistic模型)r=0.2557,xm=392.1-3-1711Anna-3-17Anna12模型检验用模型计算美国人口，与实际数据比较实际为281.4(百万)模型应用——预报美国人口加入人口数据后重新预计模型参数Logistic模型在经济领域中应用(如耐用消费品售量)阻滞增加模型(Logistic模型)r=0.2490,xm=434.0x()=306.0-3-1713Anna传染病模型问题

描述传染病传输过程

分析受感染人数改变规律

预报传染病高潮到来时刻

预防传染病蔓延伎俩

按照传输过程普通规律，用机理分析方法建立模型-3-1714Anna

已感染人数(病人)i(t)

每个病人天天有效接触(足以使人致病)人数为模型1假设若有效接触是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)建模?-3-1715Anna模型2区分已感染者(病人)和未感染者(健康人)假设1）总人数N不变，病人和健康人百分比分别为2）每个病人天天有效接触人数为,且使接触健康人致病建模~日接触率SI模型-3-1716Anna模型21/2tmii010ttm~传染病高潮到来时刻(日接触率)tmLogistic模型病人能够治愈！?t=tm,di/dt最大-3-1717Anna模型3传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人天天治愈百分比为~日治愈率建模~日接触率1/~平均感染期~一个感染期内每个病人有效接触人数，称为接触数。-3-1718Anna模型3i0i0接触数=1~阈值感染期内有效接触感染健康者人数不超出病人数1-1/i0模型2(SI模型)怎样看作模型3(SIS模型)特例idi/dt01>10ti>11-1/i0t1di/dt<0-3-1719Anna模型4传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者百分比分别为2）病人日接触率

,日治愈率,

接触数=/建模需建立两个方程-3-1720Anna模型4SIR模型无法求出解析解在相平面上研究解性质-3-1721Anna模型4消去dtSIR模型相轨线定义域相轨线11si0D在D内作相轨线图形，进行分析-3-1722Annasi101D模型4SIR模型相轨线及其分析传染病蔓延传染病不蔓延s(t)单调减相轨线方向P1s0imP1:s0>1/i(t)先升后降至0P2:s0<1/

i(t)单调降至01/~阈值P3P4P2S0-3-1723Anna模型4SIR模型预防传染病蔓延伎俩(日接触率)卫生水平(日治愈率)医疗水平传染病不蔓延条件——s0<1/预计

降低s0提升r0

提升阈值1/

降低(=/),群体免疫-3-1724Anna模型4SIR模型被传染人数预计记被传染人数百分比x<<s0i0P1i00,s01小,s01提升阈值1/降低被传染人数百分比xs0-1/=-3-1725Anna

种群弱肉强食(食饵-捕食者模型)

种群甲靠丰富天然资源生存，种群乙靠捕食甲为生，形成食饵-捕食者系统，如食用鱼和鲨鱼，美洲兔和山猫，害虫和益虫。

模型历史背景——一次世界大战期间地中海渔业捕捞量下降(食用鱼和鲨鱼同时捕捞)，不过其中鲨鱼百分比却增加，为何？食饵（甲）数量x(t),

捕食者（乙）数量

y(t)甲独立生存增加率r乙使甲增加率减小，减小量与

y成正比乙独立生存死亡率d甲使乙死亡率减小，减小量与x成正比方程(1),(2)无解析解食饵-捕食者模型(Volterra)a~捕食者掠取食饵能力b~食饵供养捕食者能力Volterra模型平衡点及其稳定性平衡点稳定性分析P点稳定性不能用近似线性方程分析p=0,q>0P:临界状态q<0P´不稳定tx(t)y(t)020.00004.00000.100021.24063.96510.22.56493.94050.300023.97633.9269………5.10009.616216.72355.9.017316.2064………9.500018.47504.04479.600019.61363.99689.700020.83113.9587用数学软件MATLAB求微分方程数值解x~y平面上相轨线计算结果（数值，图形）x(t),y(t)是周期函数，相图(x,y)是封闭曲线观察，猜测x(t),y(t)周期约为9.6xmax65.5,xmin6,ymax20.5,ymin3.9用数值积分可算出

x(t),y(t)一周期平均值：x(t)平均值约为25,

y(t)平均值约为10。食饵-捕食者模型(Volterra)

消去dt用相轨线分析点稳定性c由初始条件确定取指数x0fmf(x)x0g(y)gmy0y0在相平面上讨论相轨线图形用相轨线分析点稳定性相轨线时无相轨线以下设y2y1xQ3Q4qy1y2x1x2pyy0xx0P0x1x2Q1Q2Q1(x1,y0),Q2(x2,y0)Q3(x,y1),Q4(x,y2)相轨线退化为P点

存在x1<x0<x2,使f(x1)=f(x2)=p存在y1<y0<y2,使g(y1)=g(y2)=q相轨线是封闭曲线族xQ3Q4f(x)xx0fm0g(y)gmy0y0相轨线P~中心相轨线是封闭曲线x(t),y(t)是周期函数(周期记T)求x(t),y(t)在一周期平均值轨线中心用相轨线分析点稳定性•T2T3T4T1PT1

T2

T3

T4x(t)“相位”领先y(t)模型解释初值相轨线方向模型解释r~食饵增加率d~捕食者死亡率b~食饵供养捕食者能力捕食者数量食饵数量Pr/ad/ba~捕食者掠取食饵能力捕食者数量与r成正比,与a成反比食饵数量与d成正比,与b成反比模型解释一次大战期间地中海渔业捕捞量下降，不过其中鲨鱼百分比却在增加，为何？rr-1,dd+1捕捞战时捕捞rr-2,dd+2,2<1•••xy食饵(鱼)降低，捕食者(鲨鱼)增加自然环境

还表明：对害虫(食饵)—益虫(捕食者)系统，使用灭两种虫杀虫剂,会使害虫增加，益虫降低。食饵-捕食者模型(Volterra)缺点与改进Volterra模型改写多数食饵—捕食者系统观察不到周期震荡,而是趋向某个平衡状态,即存在稳定平衡点加Logistic项有稳定平衡点

相轨线是封闭曲线，结构不稳定——一旦离开某一条闭轨线，就进入另一条闭轨线，不恢复原状。

自然界存在周期性平衡生态系统是结构稳定，即偏离周期轨道后，内部制约使系统恢复原状。食饵-捕食者模型(Volterra)缺点与改进r1=1,N1=20,1=0.1,w=0.2,r2=0.5,2=0.18相轨线趋向极限环结构稳定相关Matlab知识命令：常微分方程符号解函数：dsolve格式：r=dsolve(‘eq1,eq2,…’,’cond1,cond2,…’,’v’)说明：对给定常微分方程（组）eq1,eq2,…中指定符号自变量v，与给定边界条件和初始条件cond1,cond2,…，求符号解（即解析解）r；若没有指定变量v，则默认变量为t；在微分方程（组）表示式eq中，大写字母D表示对自变量（设为x）微分算子：D=s/sx,D2=d2/x2,….微分算子D后面字母则表示为因变量，即待求解未知函数.-3-1740Anna相关Matlab知识初始和边界条件由字符串表示：y(a)=b,Dy(c)=d,D2y(e)=f等，分别表示若边界条件少于方程（组）阶数，则返回结果r中会出现任意常数C1,C2,…；dsolve命令最多能够接收12个输入变量（包含方程组与定解条件个数，当然能够做到输入方程个数多于12个，只要将多个方程置于一个字符串内即可）.若没有给定输出参量，则在命令窗口显示解列表.-3-1741Anna相关Matlab知识若该命令找不到解析解，则返回一警告信息，同时返回一空sym对象.这时，用户能够用命令ode23或ode45求解方程组数值解.例1：解常微分方程：（注：【程序】『输出结果』）（1）求通解；【>>s=dsolve('Dy=a*y+b')】『s=-b/a+exp(a*t)*C1』-3-1742Anna相关Matlab知识（2）初值问题【>>dsolve('Dy=y-2*t/y','y(0)=1')】『ans=(2*t+1)^(1/2)』（3）高阶方程【>>dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0','x')】『ans=4/3*cos(x)-1/3*cos(2*x)』-3-1743Anna相关Matlab知识（4）边值问题【>>dsolve('x*D2y-3*Dy=x^2','y(1)=0','y(5)=0','x')】『ans=31/468*x^4-1/3*x^3+125/468』（5）方程组问题【>>S=dsolve('Df=f+g','Dg=-f+g','f(0)=1','g(0)=2');>>%S是一个结构>>S.f,S.g】-3-1744Anna相关Matlab知识『ans=exp(t)*(2*sin(t)+cos(t))ans=exp(t)*(2*cos(t)-sin(t))』（6）无解析解问题【>>dsolve('Dy=x+y^2','y(0)=0','x')】『ans=(3^(1/2)*AiryAi(1,-x)+Airy