《第四章微分方程模型》由会员分享，可在线阅读，更多相关《第四章微分方程模型（22页珍藏版）》请在金锄头文库上搜索。
1、第四章微分方程模型一、 微分方程模型的建立在实际问题中经常需要寻求某个变量y随另一变量t的变化规律：y=y(t)，然而常常不能直接求出。有时容易建立包含变量及导数在内的关系式，即建立变量能满足的微分方程。通过求解微分方程对所研究的问题进行解释说明。因此，微分方程建模是数学建模的重要方法，微分方程模型应用也十分广泛。建立微分方程模型时，经常会遇到一些关键词，比如速率” 增长”衰变”，边际等，常涉及到导数，再结合问题所涉及的基本规律就可以得到相应的微分方程。常用微分方程建立数学模型的方法有：(1)按规律直接列方程例1 一个较热的物体置于室温为180c的房间内，该物体最初的温度是600c, 3分钟以
2、后降到500c.想知道它的温度降到 300c需要多少时间？ 10分钟以后它的温度是多少？模型建立：根据牛顿冷却(加热)定律：将温度为T的物体放入处于常温 m的介质中时，T的变化速率正比于 T与周围介质的温度差。T(t), t 0 , T的变化速率正比于 T与周围介质的温度差,设物体在冷却过程中的温度为 即dT与T- m成正比。dt建立微分方程(4.1)(4.2)?dj7 = -k(T- m),?dt?T(Q) = 60.其中参数 k0, m=18.求得一般解为 In(T m)= kt+c,或 T=m + ce-kt,tQ,1 16代入条件，求得c=42 , k= In ,最后得3211|n叫T
3、(t)=18 + 42e3 21 ,t 0.结果：(1)该物体温度降至300Qc需要8.17分钟。(2) 10分钟以后它的温度是(2)微元分析法1T(10)= 18 + 42e3In161021=25.870c该方法的基本思想是通过分析研究对象的有关变量在一个很短时间内的变化情况，寻求一些微元之间的关系式。例2 一个高为2米的球体容器里盛了一半的水，水从它的底部小孔流出，小孔的横截面积为1平方厘米试求放空容器所需要的时间2米#模型建立：首先对孔口的流速做两条假设：(1) t时刻的流速v依赖于此刻容器内水的高度h(t)。(2) 整个放水过程无能量损失。由水力学知：水从孔口流出的流量 Q为 通过孔
4、口横截面的水的体积V对时间t的变化率”即Q = dV = 0.62S 2gh(4.3)dt其中0.62是流量系数，S是孔口横截面积(单位：平方厘米)，h(t)是水面高度(单位：厘 米)，t是时间(单位：秒)当S=1平方厘米，有(4.4)dV = 0.62. 2ghdt在微小时间间隔t，t+dt 内，水面高度h(t)降至h+dh(dh0),容器中水的体积的改变量近似为(4.5)dV = - n 2dhh(t)h+dh其中r是时刻t的水面半径，右端置负号由于dh0。记 r = ,1002 - (100- h)2 =、200h- h2比较(4.4)、(4.5)两式得微分方程如下：?0.62 . 2g
5、hdt =-冗(200h- h )dh, ?ht=0=100-积分后整理得(4.6)35t =(700000 - 1000 晤 +32) 4.65、2g令h=0,求得完全排空需要约 2小时58分。(4.7)(3)模拟近似法该方法的基本思想是在不同的假设下去模拟实际的现象， 程从数学上求解或分析解的性质， 再去同实际情况作对比, 些实际的现象。例3 (交通管理问题)在交通十字路口，都会设置红绿灯.如此模拟近似所建立的微分方观察这个模型能否模拟、近似某为了让那些正行使在交叉路口或离交叉路口太近而无法停下的车辆通过路口 ，红绿灯转换中间还要亮起一段时间的黄灯.那么,黄灯应亮多长时间才最为合理呢？分析
6、 黄灯状态持续的时间包括驾驶员的反应时间，车通过交叉路口的时间以及通过刹车距离所需的时间模型建立： 记v0是法定速度，I是交叉路口的宽度，L是典型的车身长度。则车通过路口的时间为v0现来计算刹车距离。设 W为汽车的重量，卩为摩擦系数。显然，地面对汽车的摩擦力为W,其方向与运动方向相反。汽车在停车过程中，行使的距离x与时间t的关系可由下面的微分方程表示W d2xg dt2(4.8)求得,其中g为重力加速度方程(4.8)的初始条件xt=0dxdt=Vo(4.9)刹车距离就是直到速度v = 0时汽车驶过的距离，先求解二阶微分方程(4.8),对(4.8)式从0到t积分,利用条件(4.9)得dx丄二=-
7、皿 + v0 dt在条件(4.9)下对件10)式从0到t积分，得(4.10)2 02+vt(4.11)dx件10)式中令=0,可得刹车所用时间dtt0 = $ ,从而得到x(t0F。F面计算一下黄灯状态的时间A=x(t0)+I + L+TV0件12)其中T是驾驶员的反应时间(4.13)假设T = 1s, L = 4.5m, I = 9m.另外，我们取具有代表性的(1= 0.2 .当 v0 = 45、60 以及80km/h时，黄灯时间如下表所示v0(km/h)A经验法455.27s3s656.35s4s807.28s5s经验法的结果一律比我们预测的黄灯状态短些，这使人想起，许多交叉路口红绿灯的设
8、计可能使车辆在绿灯转为红灯时正处于交叉路口。一、微分方程的数值解在高等数学中，介绍了一些特殊类型微分方程的解析解法，但是大量的微分方程由于过于复杂而往往难于求出其解析解。此时只能依靠数值解法，数值解法一般只能求得微分方程的#近似解。考虑一阶常微分方程的初值问题鱼dx?y/.Vf#在区间a, b上的解，其中f(x,y)为x,y的已知函数，y为给定的初始值，将上述问题的精确解记为y(x)。数值方法的基本思想是：在解的存在区间上取 n + 1个节点a = x0 x1 x2 ? xn = b这里差hi=xi+1 -xi, i = 0,1,n称为由xi到xi+i的步长。这些hi可以不相等，但一般取成相等
9、的，这时h =。在这些节点上采用离散化方法，(通常用数值积分、微分，泰勒n展开等)将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的 yn就是上述初值问题在节点 xn上的数值解。一般说来，不同的离散 化导致不同的方法。欧拉法是解初值问题的最简单的数值方法。对(4.14)式积分可得以下积分方程x(4.15)(4.16)y(x) = yo+ jf(t, y(t)dt当x= X1时,x1y(xj= yo + /f(t, y(t)dt要得到y(xj的值，就必须计算出(4.16)式右端的积分。但积分式中含有未知函数，无法直 接计算，只好借助于数值积分。假如
10、用矩形法进行数值积分则x f(t, y(t)dt f(X。, y(x。) (x - X。)因此有y(xj y。+ f(x,y(x)(X1 - x。)=yo +hf (xo, yo) = y1利用y及f(X1,y又可以算出y(x2)的近似值：y2 = y1 +hf (X1，yJ一般地，在任意点 Xn+1 = (n + 1)h处y(x)的近似值由下式给出yn+1 = yn +hf(Xn，yn)(4.17)其中h为步长，(4.17)式亦称为显式欧拉公式。一般而言，欧拉方法计算简便，但计算精度低，收敛速度慢。可以想象，若用梯形公式计算(4.16)式右端的积分，则可望得到较高的精度。这时jf (t, y
11、(t)dt - f (Xo,y(Xo) + f (刘，y(xj)(禺-x。)将这个结果代入(4.16)并将其中的y(X1)用y1近似代替，则得1y1 = yo + hlf(Xo,yo)+ f(X1,y1)这里得到了一个含有y1的方程式，如果能从中解出y1，用它作为y(x”的近似值，可以认为比用欧拉法得出的结果要好些。仿照求y1的方法，可以逐个地求出y2,y3,。一般地当求出yn以后，要求yn+1，则可归结为解方程：yn+1 = yn+ ：f(Xn,yn)+ f (Xn+1, yn+1)(4.18)这个方法称为梯形法则。(4.18)式亦称为梯形公式。可以证明梯形公式比欧拉公式精度高，收敛速度快。
12、然而，用梯形法则求解，需要解含有yn+1的方程式，这常常很不容易。为此,? 9?.在实际计算时，可将欧拉法与梯形法则相结合，计算公式为y + hf (Xn, yn)(4.19)k = 0,1, 2,.h=yn + 2 f(Xn, yn) + f(Xn+1, yn+1)这就是先用欧拉法由(Xn, yn)得出y(Xn+1)的初始近似值yn+1，然后用(4.19)中第二式进行迭代，反复改进这个近似值，直到yn+； - yn+1 (为所允许的误差)为止，并把y,+1取作y(Xn+1)的近似值yn+1。这个方法就叫改进欧拉方法。通常把(4.18)叫做预报校正公式，其中第一式叫预报公式，第二式叫校正公式。
13、由于(4.19 )式也是显示公式，所以采用改进欧拉方法不仅计算方便，而且精度较高，收敛速度快，是人们常用的方法之一。此外，常用的方 法还有二阶、四阶龙格库塔法和线性多步法等。三、利用MATLAB求解微分方程MATLAB对常微分方程的求解是基于一阶方程进行的，通常采用龙格一一库塔方法，所对应的 MATLAB 命令为 ode (Odinary Differential Equation 的缩写)，例如 ode23、ode45、ode23s、ode23tb、ode15s、ode113等，分别用于求解不同类型的微分方程，如刚性方程和非 刚性方程等。MATLAB求解微分方程的命令如下：t,x=solve
14、r(f,tspa n, xO,opti ons)其中solver可取如ode45,ode23等函数名，f为一阶微分方程组编写的M文件名，tspan为时间矢量，可取两种形式：(a)tspan=t 0,t f时，可计算出从t 。到tf的微分方程的解；(b) tspan= t0,t1,t2,? , tm时，可计算出这些时间点上的微分方程的解。x0为微分方程的初值，options用于设定误差限(缺省时设定相对误差10 -3,绝对误差10-6), 命令为：options=odeset ( reltol ,rt, abstol ,at )，其中 rt，at:分别为设定的相对误 差和绝对误差。输出变量x记录着微分方程的解，t包含相应的时间点。下面我们按步骤给出用MATLA求解微分方程的过程。(1) 首先将常微分方程变换成一阶微分方程组。例如对以下微分方程：y(n)= f(,y,2? ,y(n-1)若令y1 = y, y2 = y,? ,yn = y(n-1)，则可得到一阶微分方程组：?1 = y2(4.20)?2 = y3?n = f(t,y1,y2,? yn)相应地可以确定初值：x0
