1、第五章第五章 微分方程模型微分方程模型常微分方程的数值解常微分方程的数值解5.1 传染病模型传染病模型常微分方程的数值解及实验常微分方程的数值解及实验在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便满足精确度要求的便于计算的表达式于计算的表达式。高数中微分方程解法在实际中基本不会直接使用高数中微分方程解法在实际中基本不会直接使用0000121212yf(x,y) xy(x ) y(x ),(),() y ,innnyxxxxxy xy xyy对
2、常微分方程 ：，其数值解是指由初始点开始的若干离散的处，即对， 求出准确值的相应近似值。因此，研究常微分方程的数值解法十分必要因此，研究常微分方程的数值解法十分必要。001i)y(xy)f(x,y , 1, 2 , 1 , 0 , xynihxi解微分方程：可用以下离散化方法求设1、用差商代替导数、用差商代替导数 若步长h较小，则有hxyhxyxy)()()( 故有公式：1-n,0,1,2,i )(),(001xyyyxhfyyiiii此即欧拉法欧拉法(向前欧拉法向前欧拉法).对应有隐式对应有隐式欧拉法欧拉法11100(,)() iiiiyyhf xyyy x2、使用数值积分、使用数值积分对方
3、程y=f(x,y), 两边由xi到xi+1积分,并利用梯形公式梯形公式, 有11()()( ,( )iixiixy xy xf x y xdx实际应用时，与欧拉公式结合使用, 2 , 1 , 0 ),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyykiiiiikiiiii的计算。然后继续下一步，取时，当满足，对于已给的精确度）（ y y 2i111i)(1)1(1kikikiyyy故有公式)(),(),(200111xyyyxfyxfhyyiiiiii梯形方法梯形方法/*trapezoid formula*/此即改进的欧拉法改进的欧拉法111(,()(,()()2iii
4、iiif xy xf xy xxx中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii1()yx).,(),2,2(),2,2(),(),22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn3、使用泰勒公式、使用泰勒公式以此方法为基础，有龙格龙格-库塔法库塔法、线性多步法线性多步法等方法1123121312(4),6(,),(,),22(,2).nnnnnnnnhyyKKKKfxyhhKfxyKKfxhyh
5、 Kh K库塔库塔三阶方法四阶龙格-库塔公式4、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为o(hk)时(k为正整数，h为步长)，称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。t，x= solver( f, ts, x0, options）ode45 ode23 ode113ode15sode23s由待解由待解方程写方程写成的成的m-文件名文件名ts=t0，tf，t0、tf为自变量为自变量的初值和终值的初值和终值函数初函数初值条件值条件自变自变量值量值函数函数值值用于设定误差限误
6、差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset(reltol,rt,abstol,at), rt, at：分别为设定的相对误差相对误差和绝对误绝对误差.t,x=ode23(f, ts, x0) 3级级2阶龙格阶龙格-库塔公式库塔公式 t,x=ode45(f, ts, x0) 5级级4阶龙格阶龙格-库塔公式库塔公式 1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:设取步长 ，从 到 用四阶龙格-库塔方法.1)0
7、(),10(2yxyxyy2.0h0 x1x求解初值问题h=0.2;ts=0:h:1;y0=1;t, x=ode45(dfun1,ts,y0);t, x, plot(t,x)function dx=dfun1(x,y)dx=y-2*x/y; 0 1.0000 0.2000 1.1832 0.4000 1.3416 0.6000 1.4832 0.8000 1.6125 1.0000 1.7321建立m-文件dfun1.m如下输入命令3、结果如图解解: 令 y1=x，y2=y1,1、建立、建立m-文件文件dfun2.m如下：如下：function dx=dfun2(t,y)dx=y(2); (1
8、-y(1)2)*y(2)-y(1);222(1)0(0)2;(0)0d xdxxxdtdtxx例例则微分方程变为一阶微分方程组：122212112(1)(0)2,(0)0yyyyyyyyt,y=ode45(dfun2,0,20,2,0);t,yplot(t,y(:,1),r-,t,y(:,2),b.-);hold on plot(y(:,1),y(:,2),co);hold offlegend(t x, t x, x y);2、取t0=0, tf=20, 输入命令：222(1)0(0)2;(0)0d xdxxxdtdtxx122212112(1)(0)2,(0)0yyyyyyyy3、结果如图t
9、,x=ode45(dfun2,0,20,2,0)plot(t,x(:,1),r-,t,x(:,2),b.-);hold on plot(x(:,1),x(:,2),co);hold offlegend(t x,t x,x y); 描述对象特征随时间描述对象特征随时间(空间空间)的演变过程的演变过程 分析对象特征的变化规律分析对象特征的变化规律 预报对象特征的未来性态预报对象特征的未来性态 研究控制对象特征的手段研究控制对象特征的手段 根据函数及其根据函数及其变化率变化率之间的关系确定函数之间的关系确定函数微分方微分方程建模程建模 根据建模目的和问题分析作出简化假设根据建模目的和问题分析作出简化
10、假设 按照内在规律或用类比法建立微分方程按照内在规律或用类比法建立微分方程5.1 传染病模型传染病模型 描述传染病的传播过程描述传染病的传播过程 分析受感染人数的变化规律分析受感染人数的变化规律 预报传染病高潮到来的时刻预报传染病高潮到来的时刻 预防传染病蔓延的手段预防传染病蔓延的手段 本世纪初, 瘟疫常在世界上某地流行, 随着人类文明的不断进步, 很多疾病, 诸如天花、霍乱已经得到有效的控制. 然而, 即使在今天, 一些贫穷的发展中国家，仍出现传染病流行的现象,医疗卫生部门的官员与专家所关注的问题是：问题提出问题提出感染疾病的人数与哪些因素有关感染疾病的人数与哪些因素有关? 不同类型传染病的
11、传播过程有不同的特点, 故不从医学的角度对各种传染病的传播过程一一进行分析，而是按一般的传播机理建立模型 由于传染病在传播的过程涉及因素较多，在分析问题的过程中，不可能通过一次假设建立完善的数学模型.思路思路：针对结果中的不合理之处结果中的不合理之处，逐步修改假设，最终得出较好的模型。先做出最简单的假设，对得出的结果进行分析，模型假设：模型假设：（1）一人得病后，久治不愈，人在传染期内不会死亡。(2) 假设每个病人每天有效接触(足以使人致病)人数为 设已感染人数设已感染人数 (病人病人) x(t), 假设是连续可微函数假设是连续可微函数()( )( )x ttx tx tt 建模建模0(0 )
12、d xxd txxtx0( )tx tx e?举个实例举个实例x=0:0.1:10;y=exp(x);plot(x,y,b-);最初只有1个病人，1个病人一天可传染1个人exp(10) = 22026被传染的机会也减少，于是将变小。若有效接触的是病人，则不能使病人数增加若有效接触的是病人，则不能使病人数增加必须区分已感染者必须区分已感染者(病人病人)和未感染者和未感染者(健康人健康人)模型缺点模型缺点问题：问题：随着时间的推移，病人的数目将无限增加，这一点与实际情况不符模型修改的关键：模型修改的关键：的变化规律的变化规律tx原因：原因： 当不考虑传染病期间的出生、死亡和迁移时，一个地区的总人数
13、可视为常数, 在传染病流行初期， 较大，因此 应为时间t的函数。随着病人的增多，健康人数减少,sidtdi1)()(tits区分未感染者(健康人)和已感染者(病人)假设1）总人数N不变，健康人和病人 的比例分别为( ), ( )s ti t 2）每个病人每天有效接触人数为, 且使接触的健康人致病建模 ()( )( )( )N i tti tNi ts tt 0)0()1(iiiidtdi 日接触率SI 模型Susceptible Infectiveteiti1111)(0Logistic 模型方程的解：方程的解：0)0()1(iiiidtdi传染病患者比例与时间t关系传染病人数的变化率与患者比
14、率i的关系染病人数由开始到高峰并逐渐达到稳定增长速度由低增至最高后降落下来对模型作进一步分析对模型作进一步分析itdiidt感染病人占一半时传染率最大！模型模型21/2tm011ln1mtitm传染病高潮到来时刻传染病高潮到来时刻 (日接触率日接触率) tm ，推迟传染高峰的到来，推迟传染高峰的到来,即改善保健措施即改善保健措施, 提高卫生水平提高卫生水平可推迟传染病高潮到来可推迟传染病高潮到来.t=tm, (i=1/2), di/dt 最大最大病人最多的一天病人最多的一天日接触率日接触率 表示该地区的表示该地区的卫生水平卫生水平， 越小卫生水平越高。越小卫生水平越高。teiti1111)(0
15、it(1)diiidt缺点：缺点：当t时，i(t) 1，这表示所有的人最终都将成为病人，这一点与实际情况也不符原因：原因：这是由假设(1)所导致，没有考虑病人可以治愈病人可以治愈及病人病发身亡病人病发身亡的情况。思考题：思考题：考虑有病人病发身亡的情况，再对模型进行修改。teiti1111)(0传染病无免疫性传染病无免疫性病人治愈成为健康人，病人治愈成为健康人， 健康人可再次被感染健康人可再次被感染增加假设增加假设SIS 模型模型3）病人平均每天治愈总病人数的比例为 日日治愈率治愈率每天治愈的病人为Ni ；病人治愈后成为仍可被感染的健康者。健康者和病人在总人数中所占的比例分别为s(t)、i(t
16、)，则： s(t)+i(t)1(1/称为传染病的平均传染期传染病的平均传染期) ()( )N i tti t 建模建模 日接触率日接触率1/ 感染期感染期0)0()1(iiiiidtdi解析法可求解该模型方程的解( )( )Ni ts tt ( )Ni tt=0(0)ii=/ 一个感染期内一个感染期内每个病人的有效接触人数，称为每个病人的有效接触人数，称为接触数接触数利用将方程写成)11 (iidtdi1-1/ idi/dt01 1i0i0011( )()1i tii 时，的增减取决于 的大小，但极限0ti 11-1/ 1( )1()i tt 0111, ( )1ii t0111, ( )1i
17、i t(1)d iiiid tlm=-阙阙 值值i0i0011( )()1i tii 时，的增减取决于 的大小，但极限0ti 11-1/ 是因为随着传染期内被传染人数占当时健康人数的比例的增加，当时的病人数所占比例也随之上升当增大时，i()也增大，)11 (iidtdi控制有效接触（隔离的效果）将最终消灭传染病。控制有效接触（隔离的效果）将最终消灭传染病。1( )0()i tt 原因：原因：感染期内有效接触使健康人数变成的病人人数不超过不超过把病人治愈的人数。模型模型2(SI模型模型)如何看作模型如何看作模型3(SIS模型模型)的特例的特例i0i0t 1di/dt 1/ i(t)先升后降至先升
18、后降至0P2: s01/ i(t)单调降至单调降至01/ 是传染病蔓延是传染病蔓延与否的与否的阈值阈值P3P4P2S01()disiii sdt预防传染病蔓延的手段预防传染病蔓延的手段 (日接触率日接触率) 卫生水平卫生水平 (日日治愈率治愈率) 医疗水平医疗水平 传染病不蔓延的条件传染病不蔓延的条件s01/ 降低降低 s0提高提高 r0 1000ris 提高阈值提高阈值 1/ 降低降低 (= / ) , 群体免疫与预防群体免疫与预防ssss00lnln0ln1000sssis0i忽略相轨线相轨线一次传染病结束后，可估计出0,ss000ln1)()(sssissi 一个感染期内一个感染期内每个
19、病人的有效接触人数，称为每个病人的有效接触人数，称为接触数接触数被传染人数的估计被传染人数的估计0ln1000sssis记被传染人数比例记被传染人数比例ssx00)211(200sxsx0)1ln(10sxx)1(200ssxxs0i0s/1P10ss i0 0, s0 1 小小, s0 1提高阈值提高阈值1/ 降低降低被传染被传染人数比例人数比例 xs0 - 1/ = 1)2(x(1)22当当 1/ 时时,即即 *1/ 1 1126621266122(101)(102)(0)10 / 4,(0)10 / 41/ 2xxxxxxxx 求解刚性方程的命令求解刚性方程的命令: ode23s, od
20、e15s 等等 (用法相同用法相同)设某城市共有n+1人，其中一人出于某种目的编造了一个谣言。该城市具有初中以上文化程度的人占总人数的一半，这些人只有1/4相信这一谣言，而其他人约有1/3会相信。又设凡相信此谣言的人每人在单位时间内传播的平均人数正比于当时尚未听说此谣言的人数，而不相信此谣言的人不传播谣言。试建立一个反映谣传情况的微分方程模型。思考题思考题1 0.99000000000000 0.00042724073268 0.000569654310241.00000000000000 0.00042731613489 0.00056975484653function df=dfrumor
21、(t,x)d1=1/2; dxin=1/4;x1=1/2; xxin=1/3;chuan=1/100;df=zeros(2,1);df=chuan*d1*dxin*(1-x(1)-x(2);chuan*x1*xxin*(1-x(1)-x(2);1(1)1(1)(0)0,(0)0ddchuanddxinxddtdxchuanxxxinxddtdx凡相信此谣言的人每人在单位时间内传播的平均人数正比: chuan初中以上文化程度的人比例:d1; 这些人相信这一谣言的比例:dxin;初中以下文化程度的人比例:x1; 这些人相信这一谣言的比例:xxin;ts=0:10:1000;x0=0,0;t,x=ode45(dfrumor,ts,x0);format longt,x
