1、微分方程模型洛阳理工学院数理部第1页1微分方程模型人口增加预测传染病模型种群模型第2页2动态动态模型模型 描述对象特征随时间描述对象特征随时间(空间空间)演变过程演变过程.分析对象特征改变规律分析对象特征改变规律.预报对象特征未来性态预报对象特征未来性态.研究控制对象特征伎俩研究控制对象特征伎俩.依据函数及其改变率之间关系确定函数依据函数及其改变率之间关系确定函数.微分微分方程方程建模建模 依据建模目标和问题分析作出简化假设依据建模目标和问题分析作出简化假设.按照内在规律或用类比法建立微分方程按照内在规律或用类比法建立微分方程.第3页3对微分方程研究方法解在很广泛条件下存在，但能用有限解析式表
2、示者极少.另辟它径：1、求数值解（近似解）；2、定性方法分析.第4页4背景背景 年年 1625 1830 1930 1960 1974 1987 1999人口人口(亿亿)5 10 20 30 40 50 60世界人口增加概况世界人口增加概况中国人口增加概况中国人口增加概况 年年 1901 1929 1953 1965 1982 1990 1995 人口人口(亿亿)4.26 5.48 6.02 7.25 10.32 11.30 12.00 12.95研究人口改变规律研究人口改变规律控制人口过快增加控制人口过快增加人口增加预测人口增加预测第5页5指数增加模型指数增加模型马尔萨斯提出马尔萨斯提出 (
3、1798)惯用计算公式惯用计算公式x(t)时刻时刻t人口人口基本假设基本假设:人口人口(相对相对)增加率增加率 r 是常数是常数今年人口今年人口 x0,年增加率年增加率 rk年后人口年后人口伴随时间增加，人口按指数规律无限增加伴随时间增加，人口按指数规律无限增加第6页6年年17901800181018201830184018501860人口3.95.37.29.612.917.123.231.4年18701880189019001910192019301940人口38.650.262.976.092.0106.5123.2131.7年19501960197019801990人口150.7179
4、.3204.0226.5251.4281.4-3-177Anna美国人口统计数据第7页指数增加模型应用及不足指数增加模型应用及不足 与与19世纪以前欧洲一些地域人口统计数据吻合世纪以前欧洲一些地域人口统计数据吻合 适合用于适合用于19世纪后迁往加拿大欧洲移民后代世纪后迁往加拿大欧洲移民后代 可用于短期人口增加预测可用于短期人口增加预测 不符合不符合19世纪后多数地域人口增加规律世纪后多数地域人口增加规律 不能预测较长久人口增加过程不能预测较长久人口增加过程1919世纪后人口数据世纪后人口数据人口增加率人口增加率r r不是常数不是常数(逐步下降逐步下降)-3-178Anna第8页阻滞增加模型阻滞
5、增加模型(Logistic模型模型)人口增加到一定数量后，增加率下降原因：人口增加到一定数量后，增加率下降原因：资源、环境等原因对人口增加阻滞作用资源、环境等原因对人口增加阻滞作用且阻滞作用随人口数量增加而变大且阻滞作用随人口数量增加而变大假设假设r固有增加率固有增加率(x很小时很小时)xm人口容量（资源、环境能容纳最大数量）人口容量（资源、环境能容纳最大数量）r是是x减函数减函数-3-179Anna第9页dx/dtx0 xmxm/2xmtx0 x(t)S形曲线形曲线,x增加先快后慢增加先快后慢x0 xm/2阻滞增加模型阻滞增加模型(Logistic模型模型)-3-1710Anna第10页参数
6、预计参数预计用指数增加模型或阻滞增加模型作人口用指数增加模型或阻滞增加模型作人口预报，必须先预计模型参数预报，必须先预计模型参数 r 或或 r,xm 利用统计数据用最小二乘法作拟合利用统计数据用最小二乘法作拟合例：美国人口数据（单位例：美国人口数据（单位百万）百万）1860 1870 1880 1960 1970 1980 1990 31.4 38.6 50.2 179.3 204.0 226.5 251.4阻滞增加模型阻滞增加模型(Logistic模型模型)r=0.2557,xm=392.1-3-1711Anna第11页-3-17Anna12第12页模型检验模型检验用模型计算年美国人口，与实
7、际数据比较用模型计算年美国人口，与实际数据比较实际为实际为281.4(百万百万)模型应用模型应用预报美国年人口预报美国年人口加入年人口数据后重新预计模型参数加入年人口数据后重新预计模型参数Logistic 模型在经济领域中应用模型在经济领域中应用(如耐用消费品售量如耐用消费品售量)阻滞增加模型阻滞增加模型(Logistic模型模型)r=0.2490,xm=434.0 x()=306.0-3-1713Anna第13页传染病模型传染病模型问题问题 描述传染病传输过程描述传染病传输过程 分析受感染人数改变规律分析受感染人数改变规律 预报传染病高潮到来时刻预报传染病高潮到来时刻 预防传染病蔓延伎俩预防
8、传染病蔓延伎俩 按照传输过程普通规律，用按照传输过程普通规律，用机理分析方法建立模型机理分析方法建立模型-3-1714Anna第14页 已感染人数已感染人数(病人病人)i(t)每个病人天天有效接触每个病人天天有效接触(足以使人致病足以使人致病)人数为人数为 模型模型1 1假设假设若有效接触是病人，则若有效接触是病人，则不能使病人数增加不能使病人数增加必须区分已感染者必须区分已感染者(病病人人)和未感染者和未感染者(健康人健康人)建模建模?-3-1715Anna第15页模型模型2 2区分已感染者区分已感染者(病人病人)和未感染者和未感染者(健康人健康人)假设假设1）总人数）总人数N不变，病人和健
9、康不变，病人和健康 人人 百分比分别为百分比分别为 2）每个病人天天有效接触人数）每个病人天天有效接触人数为为,且且使接触健康人致病使接触健康人致病建模建模 日日接触率接触率SI 模型模型-3-1716Anna第16页模型模型21/2tmii010ttm传染病高潮到来时刻传染病高潮到来时刻 (日接触率日接触率)tm Logistic 模型病人能够治愈！病人能够治愈！?t=tm,di/dt 最大最大-3-1717Anna第17页模型模型3传染病无免疫性传染病无免疫性病人治愈成病人治愈成为健康人，健康人可再次被感染为健康人，健康人可再次被感染增加假设增加假设SIS 模型模型3）病人天天治愈百分比为
10、）病人天天治愈百分比为 日日治愈率治愈率建模建模 日接触率日接触率1/平均感染期平均感染期 一个感染期内一个感染期内每个病人有每个病人有效接触人数，称为效接触人数，称为接触数接触数。-3-1718Anna第18页模型模型3i0i0接触数接触数 =1 阈值阈值感染期内感染期内有效接触感染健有效接触感染健康者人数不超出病人数康者人数不超出病人数1-1/i0模型模型2(SI模型模型)怎样看作模型怎样看作模型3(SIS模型模型)特例特例idi/dt01 10ti 11-1/i0t 1di/dt 1/i(t)先升后降至先升后降至0P2:s01/i(t)单调降至单调降至01/阈值阈值P3P4P2S0-3-
11、1723Anna第23页模型模型4SIR模型模型预防传染病蔓延伎俩预防传染病蔓延伎俩 (日接触率日接触率)卫生水平卫生水平 (日日治愈率治愈率)医疗水平医疗水平 传染病不蔓延条件传染病不蔓延条件s01/预计预计 降低降低 s0提升提升 r0 提升阈值提升阈值 1/降低降低 (=/),群体免疫群体免疫-3-1724Anna第24页模型模型4SIR模型模型被传染人数预计被传染人数预计记被传染人数百分记被传染人数百分比比x 0P:临界状态临界状态 q 0P 不稳定不稳定 第28页tx(t)y(t)020.00004.00000.100021.24063.96510.22.56493.94050.30
12、0023.97633.92695.10009.616216.72355.9.017316.20649.500018.47504.04479.600019.61363.99689.700020.83113.9587用数学软件用数学软件MATLAB求求微分方程数值解微分方程数值解xy 平面上相轨线平面上相轨线第29页计算结果（数值，图形）计算结果（数值，图形）x(t),y(t)是周期函数，相图是周期函数，相图(x,y)是封闭曲是封闭曲线线观察，猜测观察，猜测x(t),y(t)周期约为周期约为9.6xmax 65.5,xmin 6,ymax 20.5,ymin 3.9用数值积分可算出用数值积分可算出
13、 x(t),y(t)一周期平均值：一周期平均值：x(t)平均值约为平均值约为25,y(t)平均值约为平均值约为10。食饵食饵-捕食者模型捕食者模型(Volterra)第30页 消去消去dt用相轨线分析用相轨线分析 点稳定性点稳定性c 由初始条件确定由初始条件确定取指数取指数第31页x0fmf(x)x0g(y)gmy0y0在相平面上讨论相轨线图形在相平面上讨论相轨线图形用相轨线分析用相轨线分析 点稳定性点稳定性相轨线相轨线时无相轨线时无相轨线以下设以下设第32页y2y1xQ3Q4qy1y2x1x2pyy0 xx0P0 x1x2Q1Q2Q1(x1,y0),Q2(x2,y0)Q3(x,y1),Q4(
14、x,y2)相轨线相轨线退化为退化为P点点 存在x1x0 x2,使f(x1)=f(x2)=p存在y1y0y2,使g(y1)=g(y2)=q相轨线是封闭曲线族相轨线是封闭曲线族xQ3Q4f(x)xx0fm0g(y)gmy0y0相轨线相轨线P中心中心第33页相轨线相轨线是封闭曲线是封闭曲线x(t),y(t)是周期函数是周期函数(周期记周期记 T)求求x(t),y(t)在一周期平均值在一周期平均值轨线轨线中心中心用相轨线分析用相轨线分析 点稳定性点稳定性第34页T2T3T4T1PT1 T2 T3 T4x(t)“相位相位”领先领先 y(t)模型解释模型解释初值初值相轨线方向相轨线方向第35页模型解释模型
15、解释r 食饵增加率食饵增加率d 捕食者死亡率捕食者死亡率b 食饵供养捕食者能力食饵供养捕食者能力捕食者捕食者 数量数量食饵食饵数量数量Pr/ad/ba 捕食者掠取食饵能力捕食者掠取食饵能力捕食者数量与捕食者数量与r成正比成正比,与与a成反比成反比食饵食饵数量与数量与d成正比成正比,与与b成反比成反比第36页模型模型解释解释一次大战期间地中海渔业捕捞量下降，不一次大战期间地中海渔业捕捞量下降，不过其中过其中鲨鱼百分比却在增加，为何？鲨鱼百分比却在增加，为何？rr-1,dd+1捕捞捕捞战时战时捕捞捕捞rr-2,dd+2,2 s=dsolve(Dy=a*y+b)】s=-b/a+exp(a*t)*C1
16、-3-1742Anna第42页相关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第43页相关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(Df=f+g,Dg=-f+g,f(0)=1,g(0)=2);%S是
17、一个结构 S.f,S.g】-3-1744Anna第44页相关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)-3-1745Anna第45页相关Matlab知识最终一个方程无解析解，系统将其转化到两个特殊函数AiryAi，实际上没有得到解.另外，要注意这里指定x为自由变量是必要，若使用【dsolve(Dy=x+y2,y(
18、0)=0)】ans=tan(t*x(1/2)*x(1/2)把自变量误认为t，而x作为参数，所以并不是所求问题解.-3-1746Anna第46页惯用数值解指令ode45：四、五阶Runge-kutta法ode23：二、三阶Runge-kutta法t,y=ode45(odefun,tspan,y0)例：解微分方程解：执行指令：odefun=inline(y-2*t/y,t,y);t,y=ode45(odefun,0,4,1);plot(t,y,o-)-3-17Anna47第47页-3-17Anna48第48页-3-17Anna49xdata=1790:10:1990;ydata=3.9000,5.3000,7.,9.6000,12.9000,17.1000,23.,31.4000,38.6000,50.,62.9000 76.0000,92.0000,106.5000,123.,131.7000,150.7000,179.3000,204.0000,226.5000,251.4000;fun=inline(a(1)./(1+(1+(-1+a(1)./3.9)*exp(-a(2)*xdata),a,xdata);a=nlinfit(xdata,ydata,0.01,-21)第49页
