4建模与仿真 微分方程模型

很难找到系统有关变量间的直接关系容易找到变量与他们的微分（变化率）之间的关系第五章微分方程模型2/392020/2/6动态模型•描述对象特征随时间(空间)的演变过程•分析对象特征的变化规律•预报对象特征的未来性态•研究控制对象特征的手段•根据函数及其变化率之间的关系确定函数微分方程建模•根据建模目的和问题分析作出简化假设•按照内在规律或用类比法建立微分方程3/392020/2/6微分方程建模的方法和步骤建立微分方程建立函数、自变量及函数对自变量的导数之间的一种平衡关系方程求解和结果分析牛顿运动定律•应用已知的模型建立新的微分方程模型•模拟近似法•根据规律列方程•微元分析法利用已知的定律，通过取极限的方法或通过任意区域上取积分的方法得到微分方程根据实际资料或大量实验数据，提出假设，找出实际现象满足的规律，建立模型，通过模型解与实际情况的对比判断该方程能否模拟实际情况数值解法、图解法4/392020/2/6微分方程建模综合实例一问题•描述传染病的传播过程•分析受感染人数的变化规律•预报传染病高潮到来的时刻•预防传染病蔓延的手段•按照传播过程的一般规律，用机理分析方法建立模型传染病模型5/392020/2/6已感染人数(病人)i(t)•每个病人每天有效接触(足以使人致病)人数为模型1假设ttititti)()()(若有效接触的是病人，则不能使病人数增加必须区分已感染者(病人)和未感染者(健康人)建模0)0(iiidtdiitteiti0)(?6/392020/2/6sidtdi1)()(tits模型2区分已感染者(病人)和未感染者(健康人)假设1）总人数N不变，病人和健康人的比例分别为)(),(tsti2）每个病人每天有效接触人数为,且使接触的健康人致病建模ttNitstittiN)()]([)]()([0)0()1(iiiidtdi~日接触率SI模型7/392020/2/6teiti1111)(00)0()1(iiiidtdi模型21/2tmii010t11ln01itmtm~传染病高潮到来时刻(日接触率)tm1itLogistic模型病人可以治愈！?t=tm,di/dt最大8/392020/2/6模型3传染病无免疫性——病人治愈成为健康人，健康人可再次被感染增加假设SIS模型3）病人每天治愈的比例为~日治愈率ttNittitNstittiN)()()()]()([建模/~日接触率1/~感染期~一个感染期内每个病人的有效接触人数，称为接触数。0)0()1(iiiiidtdi9/392020/2/61,01,11)(i)]11([iidtdi模型3i0i0接触数=1~阈值/1)(ti形曲线增长按Sti)(感染期内有效接触感染的健康者人数不超过病人数小01i1-1/i0iiidtdi)1(模型2(SI模型)如何看作模型3(SIS模型)的特例idi/dt0110ti11-1/i0t1di/dt010/392020/2/6模型4传染病有免疫性——病人治愈后即移出感染系统，称移出者SIR模型假设1）总人数N不变，病人、健康人和移出者的比例分别为)(),(),(trtsti2）病人的日接触率,日治愈率,接触数=/建模1)()()(trtits需建立的两个方程)(),(),(trtsti11/392020/2/6ttNittitNstittiN)()()()]()([模型4SIR模型很小）通常000)0((1rrsi无法求出的解析解)(),(tsti在相平面上研究解的性质is~ttitNststtsN)()()]()([00)0(,)0(ssiisidtdsisidtdi12/392020/2/6functiony=ill(t,x)a=1;b=0.3;y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]';ts=0:50;x0=[0.02,0.98];[t,x]=ode45('ill',ts,x0);[t,x]plot(t,x(:,1),t,x(:,2)),gridplot(x(:,2),x(:,1)),grid13/392020/2/6i(t)、s(t)图形i~s图形（相轨线）14/392020/2/60011iisdsdiss000ln1)()(sssissi模型400)0(,)0(ssiisidtdsisidtdi/消去dtSIR模型}1,0,0),{(isisisD相轨线的定义域)(si相轨线11si0D在D内作相轨线的图形，进行分析)(si15/392020/2/6si101D模型4SIR模型相轨线及其分析)(si00)0(,)0(ssiisidtdsisidtdi0011iisdsdiss000ln1)()(sssissi0ln1000sssiss满足miis,/1传染病蔓延传染病不蔓延s(t)单调减相轨线的方向0,itP1s0/1imsP1:s01/i(t)先升后降至0P2:s01/i(t)单调降至01/~阈值P3P4P2S016/392020/2/6ssss00lnln模型4SIR模型预防传染病蔓延的手段(日接触率)卫生水平(日治愈率)医疗水平传染病不蔓延的条件——s01/的估计0ln1000sssis0i忽略•降低s0提高r01000ris•提高阈值1/降低(=/),群体免疫17/392020/2/6模型4SIR模型被传染人数的估计0ln1000sssis记被传染人数比例ssx00)211(200sxsx0)1ln(10sxx)1(200ssx2xxs0i0s/1P10ssi00,s01小,s01提高阈值1/降低被传染人数比例xs0-1/=18/392020/2/65.4药物在体内的分布与排除•药物进入机体形成血药浓度(单位体积血液的药物量)•血药浓度需保持在一定范围内——给药方案设计•药物在体内吸收、分布和排除过程——药物动力学•建立房室模型——药物动力学的基本步骤•房室——机体的一部分，药物在一个房室内均匀分布(血药浓度为常数)，在房室间按一定规律转移•本节讨论二室模型——中心室(心、肺、肾等)和周边室(四肢、肌肉等)19/392020/2/6中心室周边室给药排除)(0tf111)(),(Vtxtc222)(),(Vtxtc12k21k13k)()(02211131121tfxkxkxktx模型假设•中心室(1)和周边室(2),容积不变•药物在房室间转移速率及向体外排除速率，与该室血药浓度成正比•药物从体外进入中心室，在二室间相互转移,从中心室排出体外模型建立2,1~~)(~)(iVtctxiii容积浓度药量给药速率~0f2211122)(xkxktx20/392020/2/6tttteBeAtceBeAtc222111)()(1321132112kkkkk2211122121022112113121)()()()(ckckVVtcVtfckVVckktc2,1),()(itcVtxiii线性常系数非齐次方程对应齐次方程通解模型建立21/392020/2/6)()()(])()[()()(212022121101tttteeVkDtcekekVDtc0)0(,)0(,0)(21010cVDctf几种常见的给药方式1.快速静脉注射t=0瞬时注射剂量D0的药物进入中心室,血药浓度立即为D0/V12211122121022112113121)()()()(ckckVVtcVtfckVVckktc1321132112kkkkk给药速率f0(t)和初始条件22/392020/2/612211312121221131212213210122221130111)(,)(0,)(0,)(BVkkkVBAVkkkVATtVkkkkeBeAtcTtVkkeBeAtctttt0)0(,0)0(,)(2100ccktf2.恒速静脉滴注2211122121022112113121)()()()(ckckVVtcVtfckVVckktctT,c1(t)和c2(t)按指数规律趋于零药物以速率k0进入中心室0Tt23/392020/2/60010xkf)(0tx吸收室中心室000010)0()(DxxktxtkttEeBeAetc01)(1tkeDtx0100)(tkekDtxktf010100010)()(3.口服或肌肉注射相当于药物(剂量D0)先进入吸收室，吸收后进入中心室吸收室药量x0(t)2211122121022112113121)()()()(ckckVVtcVtfckVVckktcEBAcc,,0)0(,0)0(2124/392020/2/6ttBeAetctc)()(~11参数估计各种给药方式下的c1(t),c2(t)取决于参数k12,k21,k13,V1,V2t=0快速静脉注射D0,在ti(i=1,2,n)测得c1(ti)])()[()()(2121101ttekekVDtc充分大设t,由较大的用最小二乘法定A,)(,1iitct由较小的用最小二乘法定B,)(~,1iitctttAeeVkDtc)()()(1210125/392020/2/6211312kkkBAVDc101)0(011130)(dttcVkD0,,21cct1321132112kkkkkBAVkD1130ABBAk)(131321kk参数估计进入中心室的药物全部排除26/392020/2/6背景年1625183019301960197419871999人口(亿)5102030405060世界人口增长概况中国人口增长概况年19081933195319641982199019952000人口(亿)3.04.76.07.210.311.312.013.0研究人口变化规律控制人口过快增长5.6如何预报人口的增长27/392020/2/6指数增长模型——马尔萨斯提出(1798)常用的计算公式kkrxx)1(0x(t)~时刻t的人口基本假设:人口(相对)增长率r是常数trtxtxttx)()()(今年人口x0,年增长率rk年后人口0)0(,xxrxdtdxrtextx0)(trextx)()(0trx)1(0随着时间增加，人口按指数规律无限增长28/392020/2/6指数增长模型的应用及局限性•与19世纪以前欧洲一些地区人口统计数据吻合•适用于19世纪后迁往加拿大的欧洲移民后代•可用于短期人口增长预测•不符合19世纪后多数地区人口增长规律•不能预测较长期的人口增长过程19世纪后人口数据人口增