文天数学建模———微分方程模型

1、微分方程模型微分方程模型时正华 数学系2015年8月数学建模的一般步骤数学建模的一般步骤1.了解问题的实际背景，明确建模目的，收集掌握了解问题的实际背景，明确建模目的，收集掌握必要的数据资料。必要的数据资料。2.2.通过对资料的分析计算，找出起主要作用的因素，通过对资料的分析计算，找出起主要作用的因素，经必要的精练、简化，提出若干符合客观实际的经必要的精练、简化，提出若干符合客观实际的假设。假设。3.3.在所作假设的基础上，利用适当的数学工具去刻在所作假设的基础上，利用适当的数学工具去刻画各变量之间的关系，即建立模型。画各变量之间的关系，即建立模型。4.4.模型求解（包括解方程、图解、逻辑推理

2、、定理模型求解（包括解方程、图解、逻辑推理、定理证明等）。证明等）。5.5.模型的分析与检验。模型的分析与检验。1.0常微分方程基本理论常微分方程基本理论与建模方法与建模方法a. a. 常微分方程解的存在唯一性常微分方程解的存在唯一性c. c. 常微分方程的定性方法常微分方程的定性方法-相图与稳定性分析相图与稳定性分析1.0.1常微分方程的基本问题b. b. 常微分方程的求解方法常微分方程的求解方法-初等积分法与数值方法初等积分法与数值方法d. d. 常微分方程建模常微分方程建模-两个简单实例两个简单实例定义定义：一般高阶线性常微分方程具有形式通解初值问题（初值问题（IVPIVP）与特解 00

3、( , ),( ),dxf t xdtx tx),.,;(11 nnnndtxddtdxxtfdtxd),.,1 , 0(niydtxdiii 令令 ),.,;(.,11012110nnyyytfdtdyydtdyydtdya. a. 常微分方程基本概念及解的存在唯一性常微分方程基本概念及解的存在唯一性定理：定理：对初值问题（IVP）00( , )(1)( )(2)dxf t xdtx tx001212120( , )( , ):|,|.0,|( ,)( ,)|,( ,),( ,).(1)|(2)( )min( ,),max |( )|.Rf t xRt xtta xxbxLipschitzL

4、f t xf t xL xxt xt xRtthxtbhaMf xm设函数在矩形域上连续,且关于 满足条件 即存在常数使得则方程在区间上存在满足的解且解唯一.其中常微分方程解的存在唯一性常微分方程解的存在唯一性b. b. 微分方程的求解方法微分方程的求解方法-解析解与数值解解析解与数值解一阶微分方程：一阶微分方程：可分离变量可分离变量, , 齐次方程齐次方程, , 线性方程等线性方程等. .二阶齐次微分方程二阶齐次微分方程: :0axbxcx121121221212212240,40,240,(cossin)4,22rtr trtrttbacxc ec ebbacxc ec terrabacx

5、ectctbacbaa 221244,22bbacbbacrraa （1 1）解析解法）解析解法-初等积分法初等积分法解特征方程解特征方程用用MatlabMatlab求微分方程的解析解求微分方程的解析解命令：命令：dsolve(方程方程1 1, , 方程方程2 2, ,方程方程n n, , 初始条件初始条件, , 自变量自变量)记号: 在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分. 任何D后所跟的字母为因变量，自变量可以指定或由系统规则选定为确省.例如，微分方程有 应表达为：D2y=0.220d ydx用用MatlabMatlab求微分方程的解析解求微分方程的解析解例例1

6、1：求微分方程求微分方程 的通解。的通解。解：输入命令:dsolve(Du=1+u2,t)21duudt 输出结果:u=tan(t+c1)例例 2 求微分方程的特解. 15)0( , 0)0(029422yyydxdydxyd解解: : 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)输出结果:y=3e-2xsin（5x） 例例 3 求微分方程组的通解. zyxdtdzzyxdtdyzyxdtdx244354332解解: : 输入命令 ：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2

7、*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)e2tb. b. 微分方程的求解方法微分方程的求解方法-解析解与数值解解析解与数值解数学软件：数学软件：matlab, mathematica, mapleMatlab：ode45（非刚性方程非刚性方程） ode23（刚性方程刚性方程）（2 2）数值解法数值解法Mathematica：ode45（非刚

8、性方程非刚性方程） ode23（刚性方程刚性方程）用用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，x=solver（f,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和

9、绝对误差.例1. 求解下列初值问题：2,(0)1.xxx 解： （1）创建函数文件：xprim1.mfunctiony=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则微分方程变为一阶微分方程组：0)0(, 2)0()1 (1000211221221yyyyyyyy1、建立

10、m-文件vdp1000.m如下： function dy=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、结果如图050010001500200025003000-2.5-2-1.5-1-0.500.511.52例3.求解下列初值问题：11121221220.10.01 ,(0)30,0.020.04 ,(0)20.xxx xt xxxx xt x 解：（1）创建函

11、数文件：xprim2.mfunctiony=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如下： function dy=rigid(t,y) dy=

12、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,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)3、结果如图024681012-1-0.8-0.6-0.4-0.200.20.40.60.81图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.C.C. 常微分方程的定性方法常微分方程的定性方法-相图与稳定性分析相图与稳定性分析如果如果0)(limxtxt则称平衡点则称

13、平衡点x0是是稳定稳定的的.d( )(1)dxf xt称代数方程称代数方程f (x)=0的实根的实根x =x0为方程的为方程的平衡点平衡点(或奇点或奇点).它也是方程的解它也是方程的解.定义：定义：设设定性分析方法（几何方法）：定性分析方法（几何方法）：不求方程的解，而直不求方程的解，而直接分析解的（几何）性质。接分析解的（几何）性质。稳定性判别方法稳定性判别方法由于由于),)()(00 xxxfxf在讨论方程在讨论方程(1)(1)的的00d()()(2)dxfxxxt来代替来代替. .稳定性时，可用其线性近似稳定性时，可用其线性近似 易知易知x0也是方程也是方程(2)的平衡点的平衡点.(2)

14、的通解为的通解为,e)(0)(0 xCtxtxf关于关于x0是否稳定有以下结论：是否稳定有以下结论：若若, 0)(0 xf则则x0是稳定的；是稳定的；若若则则x0是不稳定的是不稳定的., 0)(0 xf 常微分方程组的平衡点及其稳定性常微分方程组的平衡点及其稳定性d( , ),d(3)d( , ).dxf x ytyg x yt代数方程组代数方程组.0),(, 0),(yxgyxf的实根的实根x =x0, y =y0称为方程称为方程(3)的的平衡点平衡点,记作记作P0(x0,y0).它也是方程它也是方程(3)的解的解.设设定义：定义：如果如果,)(lim,)(lim00ytyxtxtt则称平衡

15、点则称平衡点P0是是稳定稳定的的. 下面给出判别平衡点下面给出判别平衡点P P0 0是否稳定的是否稳定的判别准则判别准则. . 定理：定理：设设 则则（1）当）当p0且且q0时时,平衡点平衡点P0是稳定的；是稳定的；（2）当）当p0或或q0时时,平衡点平衡点P0是不稳定的是不稳定的.00()() , fgpPPxy0000()()det,()()ffPPxyqgfPPxyd.d.微分方程建模微分方程建模两个简单实例两个简单实例 在实际问题中，如果直接导出变量之间的函数关在实际问题中，如果直接导出变量之间的函数关系较为困难，但导出包含未知函数的导数或微分的关系较为困难，但导出包含未知函数的导数或

16、微分的关系式较为容易，可用建立微分方程模型。对于系式较为容易，可用建立微分方程模型。对于连续变连续变量量的问题，微分方程是十分常用的数学工具之一。的问题，微分方程是十分常用的数学工具之一。微分方程建模微分方程建模 根据函数及其根据函数及其变化率变化率之间的关系确定函数之间的关系确定函数 根据建模目的和问题分析作出根据建模目的和问题分析作出简化假设简化假设 按照内在规律（按照内在规律（机理分析机理分析）或用）或用类比法类比法建立建立微微分方程分方程例例1 1（理想单摆运动）建立理想单摆运动满足的微（理想单摆运动）建立理想单摆运动满足的微分方程，并得出分方程，并得出理想单摆理想单摆运动的周期公式。

17、运动的周期公式。 从图从图3-1中不难看出，小球所受的合力为中不难看出，小球所受的合力为mgsin，根据根据牛顿第二定律牛顿第二定律可得：可得：sinmlmg 从而得出两阶微分方程：从而得出两阶微分方程： 0sin0(0)0, (0)gl（5.0.1）这是理想单摆应这是理想单摆应满足的运动方程满足的运动方程（5.0.15.0.1）是一个两阶非线性方程，）是一个两阶非线性方程，不易求解。当不易求解。当很小时，很小时，sin，此时，此时，可考察（可考察（5.0.15.0.1）的近似线性方程：）的近似线性方程： 由此即可得出由此即可得出2gTl00(0)0, (0)gl（5.0.2） （5.0.25

18、.0.2）的解为）的解为: : (t)= 0cost gl其中其中当当时时,(t)=04Tt 42g Tl故有故有MQPmgl图图5-0.1（5.0.15.0.1）的近似方的近似方程程例例2 一根长度为一根长度为l l的金属杆被水平地夹在两端垂直的支架上，一的金属杆被水平地夹在两端垂直的支架上，一端的温度恒为端的温度恒为T T1 1，另一端温度恒为，另一端温度恒为T T2 2，（，（T T1 1、T T2 2为常数，为常数，T T1 1 T T2 2）。）。金属杆横截面积为金属杆横截面积为A A，截面的边界长度为，截面的边界长度为B B，它完全暴露在空气中，它完全暴露在空气中，空气温度为空气温

19、度为T T3 3，（，（T T3 3 10ti 11-1/ i0t 1di/dt 0模型模型43）假设传染病有免疫性假设传染病有免疫性病人治愈后即病人治愈后即移出感染系统，称为移出感染系统，称为移出者移出者. . 病人、健康人病人、健康人和移出者的比例分别为和移出者的比例分别为SIR模型模型改变假改变假设设4）假设总人数假设总人数N不变不变：1.)()()(trtsti5）病人的日接触率病人的日接触率 , , 日日治愈率治愈率 , , 接触数接触数 = / 建模建模1)()()(trtits需建立需建立的两个方程的两个方程)(),(),(trtsti).(),(),(trtstittNitti

20、tNstittiN)()()()()(模型模型4-第三次改进第三次改进SIR模型模型很小）通常000)0(1rrsi无法求出无法求出的解析解的解析解)(),(tsti在相平面在相平面 上研究解的性质上研究解的性质is ttitNststtsN)()()()(00)0(,)0(ssiisidtdsisidtdi0011iisdsdiss000ln1)()(sssissi模型模型400)0(,)0(ssiisidtdsisidtdi/消去消去dtSIR模型模型1,0,0),(isisisD相轨线相轨线 的定义域的定义域)(si相轨线相轨线11si0D在在D内作相轨线内作相轨线的图形，的图形，进行分

21、析进行分析)( sisi101D模型模型4SIR模型模型相轨线相轨线 及其分析及其分析)(si00)0(,)0(ssiisidtdsisidtdi0011iisdsdiss000ln1)()(sssissi0ln1000sssiss满足miis,/1传染病蔓延传染病蔓延传染病不蔓延传染病不蔓延s(t)单调减单调减相轨线的方向相轨线的方向0, itP1s0/1imsP1: i(t)先升后降至先升后降至0P2:i(t)单调降至单调降至0阈值阈值P3P4P2S0/10s/10s/1ssss00lnln模型模型4SIR模型模型预防传染病蔓延的手段预防传染病蔓延的手段 ( (日接触率日接触率) ) 卫生

22、水平卫生水平 ( (日日治愈率治愈率) ) 医疗水平医疗水平 传染病不蔓延的条件传染病不蔓延的条件s01/ 的估计的估计0ln1000sssis0i忽略降低降低 s0提高提高r01000ris提高阈值提高阈值 1/ 降低降低 (= / ) , 群体免疫群体免疫模型模型4SIR模型模型被传染人数的估计被传染人数的估计0ln1000sssis记被传染人数比例记被传染人数比例ssx00)211 (200sxsx0)1ln(10sxx)1(200ssx2xxs0i0s/1P10ssi0 0,s0 1 小小,s0 1提高阈值提高阈值1/ 降低降低被传染人数比例被传染人数比例 xs0- 1/ = 战争分类

23、：正规战争，游击战争，混合战争战争分类：正规战争，游击战争，混合战争只考虑双方兵力多少和战斗力强弱只考虑双方兵力多少和战斗力强弱兵力因战斗及非战斗减员而减少，因增援而增加兵力因战斗及非战斗减员而减少，因增援而增加战斗力与射击次数及命中率有关战斗力与射击次数及命中率有关建模思路和方法为用数学模型讨论社会领域的建模思路和方法为用数学模型讨论社会领域的实际问题提供了可借鉴的示例实际问题提供了可借鉴的示例第一次世界大战第一次世界大战Lanchester提出预测战役结局的模型提出预测战役结局的模型 战争模型战争模型0),(),()(0),(),()(tvyyxgtytuxyxftx一般模型一般模型每方战

24、斗减员率取决于双方的兵力和战斗力每方战斗减员率取决于双方的兵力和战斗力每方非战斗减员率与本方兵力成正比每方非战斗减员率与本方兵力成正比甲乙双方的增援率为甲乙双方的增援率为u(t),v(t)f, g取决于战争类型取决于战争类型x(t)甲方兵力，甲方兵力，y(t)乙方兵力乙方兵力模型模型假设假设模型模型)()(tvybxytuxayx正规战争模型正规战争模型甲方战斗减员率只取决于乙方的兵力和战斗力甲方战斗减员率只取决于乙方的兵力和战斗力双方均以正规部队作战双方均以正规部队作战xxprbbxg,忽略非战斗减员忽略非战斗减员假设没有增援假设没有增援00)0(,)0(yyxxbxyayxf(x,y)=

25、ay,a 乙方每个士兵的杀伤率乙方每个士兵的杀伤率a=ry py,ry 射击率，射击率， py 命中率命中率)(ty)(tx0ak0k0kbk0k正规战争模型正规战争模型为判断战争的结局，不求为判断战争的结局，不求x(t),y(t)而在而在相平面上讨论相平面上讨论x 与与y 的关系的关系00)0(,)0(yyxxbxyayxaybxdxdy2020bxaykkbxay22000yxk时平方律平方律模型模型甲方胜 0k平局0kyyxxprprabxy200乙方胜乙方胜游击战争模型游击战争模型双方都用游击部队作战双方都用游击部队作战甲方战斗减员率还随着甲方兵力的增加而增加甲方战斗减员率还随着甲方兵

26、力的增加而增加忽略非战斗减员忽略非战斗减员假设没有增援假设没有增援yrxxxxssrprddxyyxg/,),(00)0(,)0(yyxxdxyycxyxf(x,y)= cxy,c乙方每个士兵的杀伤率乙方每个士兵的杀伤率c =ry pyry射击率射击率py 命中率命中率py=sry /sxsx 甲方活动面积甲方活动面积sry 乙方射击有效面积乙方射击有效面积)(tycm0dm)(tx0m0m0m游击战争模型游击战争模型00)0(,)0(yyxxdxyycxyx00dxcymmdxcy乙方胜时000yxmyryyxrxxssrssrcdxy00线性线性律律模型模型甲方胜 0m平局 0mcddxd

27、y)(ty)(tx0乙方胜, 0n平局, 0n甲方胜, 0n00)0(,)0(yyxxbxycxyx混合战争模型混合战争模型甲方为游击部队，乙方为正规部队甲方为游击部队，乙方为正规部队020222bxcynnbxcy02002cxbxy乙方胜0n100)/(200 xy02002xsrsprxyryyxxx乙方必须乙方必须10倍于甲方的兵力倍于甲方的兵力设设x0=100,rx/ry=1/2,px=0.1,sx=1(km2),sry=1(m2)再生资源（渔业、林业等）与非再再生资源（渔业、林业等）与非再生资源（矿业等）生资源（矿业等）再生资源应适度开发再生资源应适度开发在持续稳产前提在持续稳产前

28、提下实现最大产量或最佳效益。下实现最大产量或最佳效益。问题问题及及 分析分析在在捕捞量稳定捕捞量稳定的条件下，如何控制捕的条件下，如何控制捕捞使产量最大或效益最佳。捞使产量最大或效益最佳。如果使捕捞量等于自然增长量，如果使捕捞量等于自然增长量，渔场鱼量渔场鱼量将保持不变将保持不变，则捕捞量稳定。，则捕捞量稳定。背景背景最优捕鱼问题 ExNxrxxFtx)1()()()1()()(Nxrxxftx)()()(xhxfxF记产量模型产量模型假设假设无捕捞时鱼的自然增长服从无捕捞时鱼的自然增长服从Logistic规律规律单位时间捕捞量与渔场鱼量成正比单位时间捕捞量与渔场鱼量成正比建模建模捕捞情况下渔

29、捕捞情况下渔场鱼量满足场鱼量满足不需要求解不需要求解x(t),只需知道只需知道x(t)稳定的条件稳定的条件r固有增长率固有增长率,N最大鱼量最大鱼量h(x)=Ex, E捕捞强度捕捞强度x(t)渔场鱼量渔场鱼量一阶微分方程的平衡点及其稳定性一阶微分方程的平衡点及其稳定性) 1 ()(xFx 一阶非线性（自治）方程一阶非线性（自治）方程F(x)=0的根的根x0微分方程的微分方程的平衡点平衡点000 xxxxx设设x(t)是方程的解，若从是方程的解，若从x0某邻域的任一初值出发，都有某邻域的任一初值出发，都有,)(lim0 xtxt称称x0是方程是方程(1)的的稳定平衡点稳定平衡点不求不求x(t),判断判