第1讲-微分方程模型

动态模型DynamicalModels,南京农业大学数学系2011年8月,数学建模的一般步骤,1.了解问题的实际背景，明确建模目的，收集掌握必要的数据资料。2.通过对资料的分析计算，找出起主要作用的因素，经必要的精练、简化，提出若干符合客观实际的假设。3.在所作假设的基础上，利用适当的数学工具去刻画各变量之间的关系，即建立模型。4.模型求解（包括解方程、图解、逻辑推理、定理证明等）。5.模型的分析与检验。,动态模型,描述对象特征随时间的演变过程,分析对象特征的变化规律,预报对象特征的未来性态,研究控制对象特征的手段,微分方程模型（连续）,差分方程模型（离散）,时间序列模型（随机）,常用动态模型,第一讲微分方程模型,1.0常微分方程基本理论与建模方法1.1传染病模型1.2经济增长模型1.3药物在体内的分布与排除1.4人口预测和控制1.5烟雾的扩散与消失1.6万有引力定律的发现,1.0常微分方程基本理论与建模方法,a.常微分方程解的存在唯一性,c.常微分方程的定性方法-相图与稳定性分析,1.0.1常微分方程的基本问题,b.常微分方程的求解方法-初等积分法与数值方法,d.常微分方程建模-两个简单实例,定义：一般高阶线性常微分方程具有形式通解初值问题（IVP）与特解,令,a.常微分方程基本概念及解的存在唯一性,定理：对初值问题（IVP）,常微分方程解的存在唯一性,b.微分方程的求解方法-解析解与数值解,一阶微分方程：可分离变量,齐次方程,线性方程等.,二阶齐次微分方程:,（1）解析解法-初等积分法,解特征方程,用Matlab求微分方程的解析解,命令：dsolve(方程1,方程2,方程n,初始条件,自变量),记号:在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.例如，微分方程有应表达为：D2y=0.,用Matlab求微分方程的解析解,例1：求微分方程的通解。,解：输入命令:dsolve(Du=1+u2,t),输出结果:u=tan(t+c1),解:输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x),输出结果:y=3e-2xsin（5x）,解:输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简y=simple(y)z=simple(z),输出结果：x=(c1-c2+c3+c2e-3t-c3e-3t)e2ty=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t,b.微分方程的求解方法-解析解与数值解,数学软件：matlab,mathematica,maple,Matlab：ode45（非刚性方程）ode23（刚性方程）,（2）数值解法,Mathematica：ode45（非刚性方程）ode23（刚性方程）,用Matlab软件求常微分方程的数值解,t，x=solver（f,ts,x0,options）,例1.求解下列初值问题：,解：,（1）创建函数文件：xprim1.m,functiony=xprim1(t,x)y=-x.2;,调用matlab的函数ode45求方程的数值解，并画出解曲线。,t,x=ode45(xprim1,0,1,1);plot(t,x,-,t,x,o);xlabel(Time:t0=0,tend=1);ylabel(xvalues:x(0)=1);,（2）创建求解函数文件：solplot.m（脚本文件），文件中,解:令y1=x，y2=y1,1、建立m-文件vdp1000.m如下：functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);,2、取t0=0，tf=3000，输入命令：T,Y=ode15s(vdp1000,03000,20);plot(T,Y(:,1),-),3、结果如图,例3.求解下列初值问题：,解：,（1）创建函数文件：xprim2.m,functiony=xprim2(t,x)y=x(1)-0.1*x(1)*x(2)+0.01*t;-x(2)+0.02*x(1)*x(2)+0.04*t;,调用matlab的函数ode45求方程的数值解，并画出解曲线及相图。,t,x=ode45(xprim2,0,20,30;20)plot(t,x);xlabel(Time:t0=0,tend=20);ylabel(xvalues:x1(0)=30,x2(0)=20);,（2）创建求解函数文件：solplot.m（脚本文件），文件中,解1、建立m-文件rigid.m如下：functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,2、取t0=0，tf=12，输入命令：T,Y=ode45(rigid,012,011);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3、结果如图,图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.,C.常微分方程的定性方法-相图与稳定性分析,定性分析方法（几何方法）：不求方程的解，而直接分析解的（几何）性质。,稳定性判别方法,易知x0也是方程(2)的平衡点.(2)的通解为,关于x0是否稳定有以下结论：,这个结论对于(4-1)也是成立的.,常微分方程组的平衡点及其稳定性,设,定义：,如果,则称平衡点P0是稳定的.,下面给出判别平衡点P0是否稳定的判别准则.定理：设,则（1）当p0且q0时,平衡点P0是稳定的；（2）当p0或q0时,平衡点P0是不稳定的.,d.微分方程建模两个简单实例,在实际问题中，如果直接导出变量之间的函数关系较为困难，但导出包含未知函数的导数或微分的关系式较为容易，可用建立微分方程模型。对于连续变量的问题，微分方程是十分常用的数学工具之一。,微分方程建模,根据函数及其变化率之间的关系确定函数,根据建模目的和问题分析作出简化假设,按照内在规律（机理分析）或用类比法建立微分方程,例1（理想单摆运动）建立理想单摆运动满足的微分方程，并得出理想单摆运动的周期公式。,从图3-1中不难看出，小球所受的合力为mgsin，根据牛顿第二定律可得：,这是理想单摆应满足的运动方程,（5.0.1）是一个两阶非线性方程，不易求解。当很小时，sin，此时，可考察（5.0.1）的近似线性方程：,由此即可得出,（5.0.1）的近似方程,例2一根长度为l的金属杆被水平地夹在两端垂直的支架上，一端的温度恒为T1，另一端温度恒为T2，（T1、T2为常数，T1T2）。金属杆横截面积为A，截面的边界长度为B，它完全暴露在空气中，空气温度为T3，（T3T2，T3为常数），导热系数为，试求金属杆上的温度分布T(x)，（设金属杆的导热率为）,一般情况下，在同一截面上的各点处温度也不尽相同，如果这样来考虑问题，本题要建的数学模型当为一偏微分方程。,但由题意可以看出，因金属杆较细且金属杆导热系数又较大，为简便起见，不考虑这方面的差异，而建模求单变量函数T(x)。,热传导现象机理：当温差在一定范围内时，单位时间里由温度高的一侧向温度低的一侧通过单位面积的热量与两侧的温差成正比，比例系数与介质有关。,由泰勒公式：,金属杆的微元x,x+dx在dt内由获得热量为：,同时，微元向空气散发出的热量为：,系统处于热平衡状态，故有：,所以金属杆各处温度T(x)满足的微分方程:,这是一个两阶常系数线性方程，很容易求解,1.1传染病模型,建模目的：,描述传染病的传播过程,分析受感染人数的变化规律,预报传染病高峰到来的时刻,预防传染病蔓延的手段,机理分析,建模方法：,已感染人数(病人)i(t),每个病人每天有效接触(足以使人致病)人数为,模型1,简化假设,机理分析与数学化,变量引入,数学模型,模型求解,结果分析,改进方案,若有效接触的是病人，则不能使病人数增加,必须区分已感染者(病人)和未感染者(健康人),失败原因,?,模型2第一次改进,变量与假设,1）区分已感染者(病人)和未感染者(健康人)。记病人和健康人的比例分别为。假设总人数N不变：,2）每个病人每天有效接触人数为,且使接触的健康人致病,机理分析,日接触率,SI模型,定性分析,Logistic模型,模型的解析求解,Logistic模型,令i=1/2得,结果分析,(1)(日接触率)tm,Logistic模型,病人可以治愈！,?!,失败原因：,模型3第二次改进,3）病人治愈成为健康人，病人每天治愈的比例为。,增加假设,SIS模型,4）健康人对传染病无免疫性，可再次被感染。,日治愈率,建模,日接触率,1/感染期,一个感染期内每个病人的有效接触人数，称为接触数。,模型3,接触数=1阈值,感染期内有效接触感染的健康者人数不超过病人数,模型2(SI模型)是模型3(SIS模型)的特例,模型4,3）假设传染病有免疫性病人治愈后即移出感染系统，称为移出者.病人、健康人和移出者的比例分别为,SIR模型,改变假设,4）假设总人数N不变：,5）病人的日接触率,日治愈率,接触数=/,建模,需建立的两个方程,模型4-第三次改进,SIR模型,模型4,SIR模型,相轨线的定义域,在D内作相轨线的图形，进行分析,模型4,SIR模型,相轨线及其分析,s(t)单调减相轨线的方向,P1:i(t)先升后降至0,P2:i(t)单调降至0,阈值,模型4,SIR模型,预防传染病蔓延的手段,(日接触率)卫生水平,(日治愈率)医疗水平,传染病不蔓延的条件s0T,c1(t)和c2(t)按指数规律趋于零,3.口服或肌肉注射,相当于药物(剂量D0)先进入吸收室，吸收后进入中心室,吸收室药量x0(t),参数估计,各种给药方式下的c1(t),c2(t)取决于参数k12,k21,k13,V1,V2,t=0快速静脉注射D0,在ti(i=1,2,n)测得c1(ti),由较大的用最小二乘法定A,由较小的用最小二乘法定B,参数估计,1.4人口预测和控制,年龄分布对于人口预测的重要性,只考虑自然出生与死亡，不计迁移,人口发展方程,人口发展方程,一阶偏微分方程,人口发展方程,已知函数（人口调查）,生育率（控制人口手段）,生育率的分解,总和生育率,h生育模式,人口发展方程和生育率,总和生育率控制生育的多少,生育模式控制生育的早晚和疏密,正反馈系统,滞后作用很大,人口指数,1）人口总数,2）平均年龄,3）平均寿命,t时刻出生的人，死亡率按(r,t)计算的平均存活时间,4）老龄化指数,控制生育率,控制N(t)不过大,控制(t)不过高,1.5烟雾的扩散与消失,现象和问题,炮弹在空中爆炸，烟雾向四周扩散，形成圆形不透光区域。,不透光区域不断扩大，然后区域边界逐渐明亮，区域缩小，最后烟雾消失。,建立模型描述烟雾扩散和消失过程，分析消失时间与各因素的关系。,问题分析,无穷空间由瞬时点源导致的扩散过程，用二阶偏微分方程描述烟雾浓度的变化。,观察的烟雾消失与烟雾对光线的吸收，以及仪器对明暗的灵敏程度有关。,模型假设,1）烟雾在无穷空间扩散，不受大地和风的影响；扩散服从热传导定律。,2）光线穿过烟雾时光强的减少与烟雾浓度成正比；无烟雾的大气不影响光强。,3）穿过烟雾进入仪器的光线只有明暗之分，明暗界限由仪器灵敏度决定。,模型建立,1）烟雾浓度的变化规律,热传导定律：单位时间通过单位法向面积的流量与浓度梯度成正比,曲面积分的奥氏公式,1）烟雾浓度的变化规律,初始条件,Q炮弹释放的烟雾总量,单位强度的点源函数,对任意t,C的等值面是球面x2+y2+z2=R2;RC,仅当t,对任意点(x,y,z),C0,1）烟雾浓度的变化规律,2）穿过烟雾光强的变化规律,光强的减少与烟雾浓度成正比,3）仪器灵敏度与烟雾明暗界限,烟雾浓度连续变化,烟雾中光强连续变化,设光源在z=-,仪器在z=,则观测到的明暗界限为,不透光区域边界,4）不透光区域边界的变化规律,对任意t,不透光区域边界是圆周,不透光区域边界半径,结果分析,观测到不透光区域边界达到最大的时刻t1，可以预报烟雾消失的时刻t2,1.6万有引力定律的发现,背景,航海业发展,天文观测精确,“地心说”动摇,哥白尼：“日心