数学建模-常微分方程模型及差分模型

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

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

分析对象特征的变化规律.

预报对象特征的未来性态.

研究控制对象特征的手段.

根据函数及其变化率之间的关系确定函数.微分方程建模

根据建模目的和问题分析作出简化假设.

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

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

年19011929195319651982199019952000人口(亿)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年195019601970198019902000人口150.7179.3204.0226.5251.4281.4美国人口统计数据2012-3-177精选课件指数增长模型的应用及局限性

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

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

可用于短期人口增长预测

不符合19世纪后多数地区人口增长规律

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

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

描述传染病的传播过程

分析受感染人数的变化规律

预报传染病高潮到来的时刻

预防传染病蔓延的手段

按照传播过程的一般规律，用机理分析方法建立模型2012-3-1714精选课件

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

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

,日治愈率,

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

i(t)单调降至01/~阈值P3P4P2S02012-3-1723精选课件模型4SIR模型预防传染病蔓延的手段(日接触率)卫生水平(日治愈率)医疗水平传染病不蔓延的条件——s0<1/的估计

降低s0提高r0

提高阈值1/

降低(=/),群体免疫2012-3-1724精选课件模型4SIR模型被传染人数的估计记被传染人数比例x<<s0i0P1i00,s01小,s01提高阈值1/降低被传染人数比例xs0-1/=2012-3-1725精选课件

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

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

模型的历史背景——一次世界大战期间地中海渔业的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但是其中鲨鱼的比例却增加，为什么？精选课件食饵（甲）数量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.200022.56493.94050.300023.97633.9269………5.10009.616216.72355.20009.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后面的字母则表示为因变量，即待求解的未知函数.2012-3-1740精选课件相关Matlab知识初始和边界条件由字符串表示：y(a)=b,Dy(c)=d,D2y(e)=f等，分别表示若边界条件少于方程（组）的阶数，则返回的结果r中会出现任意常数C1,C2,…；dsolve命令最多可以接受12个输入变量（包括方程组与定解条件个数，当然可以做到输入的方程个数多于12个，只要将多个方程置于一个字符串内即可）.若没有给定输出参量，则在命令窗口显示解列表.2012-3-1741精选课件相关Matlab知识若该命令找不到解析解，则返回一警告信息，同时返回一空的sym对象.这时，用户可以用命令ode23或ode45求解方程组的数值解.例1：解常微分方程：（注：【程序】『输出结果』）（1）求的通解；【>>s=dsolve('Dy=a*y+b')】『s=-b/a+exp(a*t)*C1』2012-3-1742精选课件相关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)』2012-3-1743精选课件相关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】2012-3-1744精选课件相关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)+AiryB