《数学建模-常微分方程模型及差分模型》由会员分享，可在线阅读，更多相关《数学建模-常微分方程模型及差分模型（49页珍藏版）》请在金锄头文库上搜索。
1、微分方程模型,洛阳理工学院数理部,1,微分方程模型,人口增长的预测 传染病模型 种群模型,2,动态模型,描述对象特征随时间(空间)的演变过程.,分析对象特征的变化规律.,预报对象特征的未来性态.,研究控制对象特征的手段.,根据函数及其变化率之间的关系确定函数.,微分方程建模,根据建模目的和问题分析作出简化假设.,按照内在规律或用类比法建立微分方程.,3,对微分方程的研究方法,解在很广泛的条件下存在，但能用有限解析式表达者很少. 另辟它径： 1、求数值解（近似解）； 2、定性方法分析.,4,背景,世界人口增长概况,中国人口增长概况,研究人口变化规律,控制人口过快增长,人口增长的预测,5,指数增长
2、模型马尔萨斯提出 (1798),常用的计算公式,x(t) 时刻t的人口,基本假设 : 人口(相对)增长率 r 是常数,今年人口 x0, 年增长率 r,k年后人口,随着时间增加，人口按指数规律无限增长,6,2012-3-17,7,Anna,美国人口统计数据,指数增长模型的应用及局限性,与19世纪以前欧洲一些地区人口统计数据吻合,适用于19世纪后迁往加拿大的欧洲移民后代,可用于短期人口增长预测,不符合19世纪后多数地区人口增长规律,不能预测较长期的人口增长过程,19世纪后人口数据,2012-3-17,8,Anna,阻滞增长模型(Logistic模型),人口增长到一定数量后，增长率下降的原因：,资源
3、、环境等因素对人口增长的阻滞作用,且阻滞作用随人口数量增加而变大,假设,r固有增长率(x很小时),xm人口容量（资源、环境能容纳的最大数量）,2012-3-17,9,Anna,x(t)S形曲线, x增加先快后慢,阻滞增长模型(Logistic模型),2012-3-17,10,Anna,参数估计,用指数增长模型或阻滞增长模型作人口 预报，必须先估计模型参数 r 或 r, xm,利用统计数据用最小二乘法作拟合,例：美国人口数据（单位百万）,阻滞增长模型(Logistic模型),2012-3-17,11,Anna,2012-3-17,Anna,12,模型检验,用模型计算2000年美国人口，与实际数据
4、比较,实际为281.4 (百万),模型应用预报美国2010年的人口,加入2000年人口数据后重新估计模型参数,Logistic 模型在经济领域中的应用(如耐用消费品的售量),阻滞增长模型(Logistic模型),2012-3-17,13,Anna,传染病模型,问题,描述传染病的传播过程,分析受感染人数的变化规律,预报传染病高潮到来的时刻,预防传染病蔓延的手段,按照传播过程的一般规律，用机理分析方法建立模型,2012-3-17,14,Anna,已感染人数 (病人) i(t),每个病人每天有效接触(足以使人致病)人数为,模型1,假设,若有效接触的是病人，则不能使病人数增加,建模,?,2012-3-
5、17,15,Anna,模型2,区分已感染者(病人)和未感染者(健康人),假设,1）总人数N不变，病人和健康 人的 比例分别为,2）每个病人每天有效接触人数为, 且使接触的健康人致病,建模, 日 接触率,SI 模型,2012-3-17,16,Anna,模型2,tm传染病高潮到来时刻, (日接触率) tm,病人可以治愈！,?,t=tm, di/dt 最大,2012-3-17,17,Anna,模型3,传染病无免疫性病人治愈成为健康人，健康人可再次被感染,增加假设,SIS 模型,3）病人每天治愈的比例为, 日治愈率,建模, 日接触率,1/ 平均感染期, 一个感染期内每个病人的有效接触人数，称为接触数。
6、,2012-3-17,18,Anna,模型3,接触数 =1 阈值,感染期内有效接触感染的健康者人数不超过病人数,模型2(SI模型)如何看作模型3(SIS模型)的特例,2012-3-17,19,Anna,模型4,传染病有免疫性病人治愈后即移出感染系统，称移出者,SIR模型,假设,1）总人数N不变，病人、健康人和移出者的比例分别为,2）病人的日接触率 , 日治愈率, 接触数 = / ,建模,需建立 的两个方程,2012-3-17,20,Anna,模型4,SIR模型,2012-3-17,21,Anna,模型4,SIR模型,相轨线 的定义域,在D内作相轨线 的图形，进行分析,2012-3-17,22,
7、Anna,模型4,SIR模型,相轨线 及其分析,s(t)单调减相轨线的方向,P1: s01/ i(t)先升后降至0,P2: s01/ i(t)单调降至0,1/阈值,2012-3-17,23,Anna,模型4,SIR模型,预防传染病蔓延的手段, (日接触率) 卫生水平,(日治愈率) 医疗水平,传染病不蔓延的条件s01/, 的估计,降低 s0,提高 r0,提高阈值 1/,2012-3-17,24,Anna,模型4,SIR模型,被传染人数的估计,记被传染人数比例, 小, s0 1,提高阈值1/降低被传染人数比例 x,s0 - 1/ = ,2012-3-17,25,Anna,种群的弱肉强食(食饵-捕食
8、者模型),种群甲靠丰富的天然资源生存，种群乙靠捕食甲为生，形成食饵-捕食者系统，如食用鱼和鲨鱼，美洲兔和山猫，害虫和益虫。,模型的历史背景一次世界大战期间地中海渔业的捕捞量下降(食用鱼和鲨鱼同时捕捞)，但是其中鲨鱼的比例却增加，为什么？,食饵（甲）数量 x(t), 捕食者（乙）数量 y(t),甲独立生存的增长率 r,乙使甲的增长率减小，减小量与 y成正比,乙独立生存的死亡率 d,甲使乙的死亡率减小，减小量与 x成正比,方程(1),(2) 无解析解,食饵-捕食者模型(Volterra),a 捕食者掠取食饵能力,b 食饵供养捕食者能力,Volterra模型的平衡点及其稳定性,平衡点,稳定性分析,P
9、点稳定性不能用近似线性方程分析,p =0, q 0 P: 临界状态,q 0 P 不稳定,用数学软件MATLAB求微分方程数值解,xy 平面上的相轨线,计算结果（数值，图形）,x(t), y(t)是周期函数，相图(x,y)是封闭曲线,x(t), y(t)的周期约为9.6,xmax 65.5, xmin 6, ymax 20.5, ymin 3.9,用数值积分可算出 x(t), y(t)一周期的平均值： x(t)的平均值约为25, y(t)的平均值约为10。,食饵-捕食者模型(Volterra),用相轨线分析 点稳定性,c 由初始条件确定,在相平面上讨论相轨线的图形,用相轨线分析 点稳定性,相轨线
10、,时无相轨线,以下设,相轨线,P中心,相轨线是封闭曲线,求x(t), y(t) 在一周期的平均值,轨线中心,用相轨线分析 点稳定性,x(t) 的“相位”领先 y(t),模型解释,初值,相轨线的方向,模型解释,r 食饵增长率,d 捕食者死亡率,b 食饵供养捕食者能力,捕食者 数量,食饵数量,a 捕食者掠取食饵能力,捕食者数量与r成正比, 与a成反比,食饵数量与d成正比, 与b成反比,模型解释,一次大战期间地中海渔业的捕捞量下降，但是其中鲨鱼的比例却在增加，为什么？,rr-1, dd+1,捕捞,战时捕捞,rr-2, dd+2 , 2 1,食饵(鱼)减少， 捕食者(鲨鱼)增加,自然环境,还表明：对害
11、虫(食饵)益虫(捕食者)系统，使用灭两种虫的杀虫剂, 会使害虫增加，益虫减少。,食饵-捕食者模型(Volterra)的缺点与改进,Volterra模型,多数食饵捕食者系统观察不到周期震荡,而是趋向某个平衡状态,即存在稳定平衡点,有稳定平衡点,相轨线是封闭曲线，结构不稳定一旦离开某一条闭轨线，就进入另一条闭轨线，不恢复原状。,自然界存在的周期性平衡生态系统是结构稳定的，即偏离周期轨道后，内部制约使系统恢复原状。,食饵-捕食者模型(Volterra)的缺点与改进,r1=1, N1=20, 1=0.1, w=0.2, r2=0.5, 2=0.18,相轨线趋向极限环,相关Matlab知识,命令：常微分
12、方程的符号解 函数：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-17,40,Anna,相关Matlab知识,初始和边界条件由字符串表示：y(a)=b,Dy(c)=d,D2y(e)=f等，分别表示 若边
13、界条件少于方程（组）的阶数，则返回的结果r中会出现任意常数C1,C2,； dsolve命令最多可以接受12个输入变量（包括方程组与定解条件个数，当然可以做到输入的方程个数多于12个，只要将多个方程置于一个字符串内即可）. 若没有给定输出参量，则在命令窗口显示解列表.,2012-3-17,41,Anna,相关Matlab知识,若该命令找不到解析解，则返回一警告信息，同时返回一空的sym对象. 这时，用户可以用命令ode23或ode45求解方程组的数值解. 例1：解常微分方程： （注：【程序】输出结果） （1）求 的通解； 【 s=dsolve(Dy=a*y+b)】 s = -b/a+exp(a*
14、t)*C1,2012-3-17,42,Anna,相关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-17,43,Anna,相关Matlab知识,（4）边值问题 【 dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x)】 ans = 31/468*x4-1/3*x3+125/468 （5）方程组问题 【 S=dsolve(
15、Df=f+g,Dg=-f+g,f(0)=1,g(0)=2); %S是一个结构 S.f,S.g】,2012-3-17,44,Anna,相关Matlab知识,ans = exp(t)*(2*sin(t)+cos(t) ans = exp(t)*(2*cos(t)-sin(t) （6）无解析解问题 【 dsolve(Dy=x+y2,y(0)=0,x)】 ans = (3(1/2)*AiryAi(1,-x)+AiryBi(1,-x)/(3(1/2)*AiryAi(-x)+AiryBi(-x),2012-3-17,45,Anna,相关Matlab知识,最后一个方程无解析解，系统将其转化到两个特殊函数AiryAi，事实上没有得到解. 另外，要注意这里指定x为自由变量是必要的，若使用 【 dsolve(Dy=x+y2,y(0)=0)】 ans = tan(t*x(1/2)*x(1/2) 把自变量误认为t，而x作为参数，所以并不是所求问题的解.,2012-3-17,46,Anna,常用数值解指令,ode45：四、五阶Runge-kutta法 ode23：二、三阶Runge-kutta法
