第4章微分方程

1、第4章 微分方程本章通过对一些简化的应用问题的分析，利用学习过的高等数学知识来解决这些实际问题。并且能够通过本章的学习，使读者学习和体会怎样收集和利用现有的一些数据，根据实际问题给出合理假设，利用这些数据和假设建立微分模型，然后通过MATLAB进行模拟给出实际问题的解决办法和控制,并在补充知识里介绍了这一模型里面的三个经典模型。4.1 引例 追线问题现实生活中，有许多实际问题都可以用微分方程来描述。例如我缉私舰雷达发现距c km处有一艘走私船正以匀速a km/h沿直线行使，缉私舰立即以最大的速度b km/h追赶，若用雷达进行跟踪，保持船的瞬时速度方向始终指向走私船，试求缉私舰追逐路线和追上的时

2、间。选取走私船逃跑方向为y轴方向，缉私舰在位置发现走私船在处，显然缉私舰、走私船的大小比他们运动的范围小得多，可视为质点。设在缉私舰发现走私船时算起的时间为，走私船到达点，缉私舰到（图4-1），因直线DR与路线相切，由几何关系得到 或 (4.1)为消去，先把(4.1)对求导，得到 (4.2)代入得到 (4.3)这里有负号是因为随的减小而增大，结合(4.2)和(4.3)，得到追线的微分方程 (4.4)其中，上式不显含，令及，则上式可化为 从而 (4.5)要继续求是的怎样一个函数，必须进一步确定。(1) 若，从而，积分式(4.5)得 当时，即走私船被缉私舰捕捉前所跑过的距离为，所用的时间是。(2)

3、 若，即，由(4.5)可得 显然不能取零，即缉私舰不可能追上走私船的。(3) 若，即，显然缉私舰不可能追上走私船。4.2 微分方程(组)简介我们知道微分方程就是联系着自变量、未知函数以及它的导数的关系式。微分方程可以分为常微分方程和偏微分方程两种。例如 (4.6) (4.7) (4.8)这里(4.6)是二阶常微分方程，而(4.8)是含有三个自变量的二阶偏微分方程。一般的阶常微分方程具有形式 (4.9)这里是的已知函数，而且一定含有,是未知函数,是自变量。4.2.1 线性微分方程和非线性微分方程如果方程(4.9)的两端为及的一次有理整式，则称(4.9)为阶线性微分方程。例如(4.6)是二阶线性微

4、分方程。一般阶线性微分方程具有形式 (4.10)这里都是的已知函数。不是线性方程的方程称为非线性方程，例如 (4.11)4.2.2 微分方程的解如果函数代入方程(4.9)后，能使它变为恒等式，则称函数为方程(4.9)的解。如果关系式决定的隐函数为方程(4.9)的解，我们称是方程(4.9)的隐式解。解微分方程的一般形式有变量分离方程、齐次方程、恰当方程、线性方程、伯努力方程以及常数变易法。4.3 求解微分方程(组)的MATLAB命令dsolve ('方程1', '方程2', '方程n', '初始条件', '自变量')

5、记号:在表达微分方程中，用字母D表示求微分，D2，D3等表示求高阶微分，任何D后所跟的字母为因变量，自变量可以指定或由系统规定而缺省。例1：求 的通解.>>u=dsolve('Du=1+u2','t')ans= u = tan(t+c)例2：求微分方程的特解.>> y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') ans=y =3e-2xsin（5x）例3：求微分方程组的通解. >>x,y,z=dsolve('Dx=2*x

6、-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)ans=x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz = (-c1e-4t+c2e-4t+c1-c2+c3)e2t 4.4 微分方程的数值解4.4.1常微分方程数值解的定义在生产和科研中所处理的微分方程往往

7、很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的。对常微分方程其数值解是由初始点开始的若干离散的处,即对求出准确值的相应近似值。4.4.2 建立数值解法的一些途径设可以用离散化方法求解微分方程4.4.2.1 用差商代替导数若步长较小，则有因此有公式此即欧拉法。4.4.2.2 使用数值积分对方程，两边由到积分，并利用梯形公式，有：故有公式实际应用时，和欧拉公式结合使用对于已给的精确度，当满足时,取，然后继续下一步的计算,此乃改进的欧拉法。4.4.2.3 使用

8、泰勒公式以此方法为基础,有龙格-库塔法、线性多步等方法。4.4.2.4 数值公式的精度当一个数值公式的截断误差可表示为时（为正整数，为步长），称它是一个阶公式, 越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。4.5 用Matlab软件求常微分方程的数值解t,x= ode45（'f',ts,x0,options）运用组合的4/5阶龙格-库塔-芬尔格算法t,x= ode23（'f',ts,x0,options）组合的2/3阶龙格-库塔-芬尔格算法t,x= ode11

9、3（'f',ts,x0,options）t,x= ode15s（'f'ts,x0,options）t,x= ode23s（'f',ts,x0,options）: 自变量值，: 函数值，：由待解方程写成的m-文件名，为自变量的初值和终值，为函数的初值，options: 用于设定误差限(缺省时设定相对误差10-3，绝对误差10-6)，命令为：options=odeset('reltol',rt, 'abstol',at) rt, at 分别为设定的相对误差和绝对误差。注意：1、 解个未知函数的方程组时，和均为维向量，m

10、-文件中的待解方程组应以的分量形式写成。2、 使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组。例4：解 令，则微分方程变为微分方程组1、 建立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、 取，输入命令T,Y=ode15s('vdp1000',0 3000,2 0); plot(T,Y(:,1),'-')3、 结果如图例5：解微分方程组.解1、建立m-文件rigid.m如下： fu

11、nction dy=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);3、 取，输入命令T,Y=ode45('rigid',0 12,0 1 1); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、结果如下图（图中，的图形为实线，的图形为“*”线，的图形为“+”线）4.6 例题4.6.1 导弹追踪问题：设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始

12、终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？解法一（解析法）假设导弹在时刻的位置为，乙舰位于，由于导弹头始终对准乙舰，故此时直线PQ就是导弹的轨迹曲线弧OP在点P处的切线，即有 (4.12)又根据题意，弧OP的长度为|AQ|的5倍，即 (4.13)由(4.12)，(4.13)消去整理得到模型 (4.14)初始条件为解即为导弹的运动轨迹：当时，即当乙舰航行到时被导弹击中。被导弹击中时间为，若，则在处被击中。解法二(数值解)令，将方程（4.14）化为一阶微分方程组。1. 建立m-文件eq1.m fu

13、nction dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(1-x); 2. 取x0=0，xf=0.9999，建立主程序ff6.m如下： x0=0，xf=0.9999 x,y=ode15s('eq1',x0 xf,0 0); plot(x,y(:,1),b.') hold on y=0:0.01:2; plot(1,y,b*')结论：导弹大致在（1，0.2）处击中乙舰解法三(建立参数方程求数值解)设时刻乙舰的坐标为，导弹的坐标为.1. 设导弹速度恒为，则 （4.15）2. 由于弹头

14、始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，即： ， （4.16）消去得： （4.17）3因乙舰以速度沿直线运动，设v0=1，则。因此导弹运动轨迹的参数方程为： (4.18)4. 解导弹运动轨迹的参数方程（1） 建立m-文件eq2.m如下： function dy=eq2(t,y) dy=zeros(2,1); dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2); dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2);（2） 取，建立主程序chase2.m如下： t,y=ode45('eq2',0 2,0 0);

15、Y=0:0.01:2; plot(1,Y,'-')， hold on plot(y(:,1),y(:,2),'*')5. 结果见下图4.6.2 慢跑者与狗：一个慢跑者在平面上沿椭圆以恒定的速率跑步，设椭圆方程为： ，,突然有一只狗攻击他。这只狗从原点出发，以恒定速率跑向慢跑者，狗的运动方向始终指向慢跑者。分别求出时狗的运动轨迹。1. 模型建立设时刻慢跑者的坐标为，狗的坐标为。则，狗从出发，与导弹追踪问题类似，建立狗的运动轨迹的参数方程： (4.19)2. 模型求解(1) 时，建立m-文件eq3.m如下： function dy=eq3(t,y) dy=zeros

16、(2,1); dy(1)=20*(10+20*cos(t)-y(1)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(2)=20*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);(2) 取，建立主程序chase3.m如下：t0=0;tf=10; t,y=ode45('eq3',t0 tf,0 0); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,'-') hold o

17、n plot(y(:,1),y(:,2),'*')在chase3.m，不断修改的值,分别取tf=5, 2.5, 3.5,至3.15时，狗刚好追上慢跑者。(3) 时建立m-文件eq4.m如下： function dy=eq4(t,y)dy=zeros(2,1);dy(1)=5*(10+20*cos(t)-y(1)/sqrt(10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);dy(2)=5*(20+15*sin(t)-y(2)/sqrt(10+20*cos(t)- y(1)2+(20+15*sin(t)-y(2)2);(4) 取，建立主程序chase4

18、.m如下： t0=0;tf=10;t,y=ode45('eq4',t0 tf,0 0);T=0:0.1:2*pi;X=10+20*cos(T);Y=20+15*sin(T);plot(X,Y,'-')hold onplot(y(:,1),y(:,2),'*')在chase3.m，不断修改的值，分别取tf=20, 40, 80,可以看出，狗永远追不上慢跑者。4.7 练习1、意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其

19、捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？ 年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.72、一个小孩借助长度为a的硬棒，拉或推某玩具.此小孩沿某曲线行走，计算并画出玩具的轨迹。3、讨论资金积累、国民收入、与人口增长的关系。（1）若国民平均收入x与按人口平均资金积累y成正比，说明仅当总资金积累的相对增长率k大于人口的相对增长率r时，国民平均收入才是增长的。（2）作出k(x)和r(x)的示意图，

20、分析人口激增会引起什么后果。4.8 补充知识微分方程模型是实际问题中很重要的一种模型，其中有几类很典型的微分方程模型1. 人口问题，利用数学模型描述、预测我国人口的变化规律，并提出相应的控制方案。2. 生物种群问题， 利用数学模型描述两种生物种群相互作用的规律，并解释为什么在战争期间捕鱼的总数减少而掠肉鱼的比例却会上升这一现象。3. 优捕鱼策略， 为了保护人类赖以生存的自然环境，可再生资源（如渔业、林业资源）的开发必须适度。一种合理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。若用表示时刻某范围内一种群的数量或密度，当种群数量较大时，我们将看作的连续函数，则的变化与出生、死

21、亡、迁入和迁出等因素有关。若用分别表示种群的出生率、死亡率、迁入率、迁出率，则种群数量或密度变化的一般模型是： (4.20)以下介绍的几个模型都是根据这个原理建立的。4.8.1 马尔萨斯(Malthus)模型考虑一个国家或地区的人口总数随时间变化的情况。记为时刻该国家或地区的人口总数，对一个国家而言,迁入和迁出人数相对较小，故略去迁移对人口变化的影响，即人口变化仅与出生率和死亡率有关。进一步假设每一个社会成员的死亡与生育水平相同，这样就可以设，即人口的出生率和死亡率之差与总人口成正比，这个常数称为自然增长率，它可以通过人口统计数据得到。这就得到了人口增长的马尔萨斯模型 (表4-1) 。表4-1

22、 我国人口预测模型年份马尔萨斯模型预测值/亿洛杰斯蒂克模型预测值/亿离散模型预测值/亿实际统计值/亿198210.154110.154110.154110.1541198310.297210.256410.305810.2495198410.442410.359410.456110.3475198510.589610.463110.606810.4532198610.738910.567310.761510.5721198710.890310.672110.923010.7240198811.043910.777511.096110.8978198911.119610.883511.28121

23、1.0676199011.357510.990111.475711.3368200013.064212.087113.5620201015.027413.235715.0905202017.285614.427616.7964203019.883215.652918.4940204022.871116.900919.7519205026.308118.159521.2622 (4.21)方程(4.21)的解是，当时,人口将以指数规律增长。在实际应用中，一般以年为间隔考察人口的变化情况，即取,这样就得到以后各年人口的总数为,这表明人口以公比为的等比级数的速度增长，这就是马尔萨斯提出“人口以几何级

24、数增长”的理论基础。下面用马尔萨斯模型估计我国人口变化情况。为了便于对比，我们取1982年人口普查时得到的人口总数为初始值，即亿，自然增长率，用公式估计以后各年我国总人口的变化，其结果如表4-1所示，从表4-1可以看到，在1983年到1990年的8年中，用马尔萨斯模型预测的相对误差均在2%以下，这表明模型比较准确的预测了短期内人口变化的规律。马尔萨斯中的关键假设是自然增长率仅与人口出生率和死亡率有关，且是常数。这一假设使该模型简单实用，但这一假设也导致了人口无限制的增长，显然用该模型来作为长期人口预测是不合理的，需要进行修改。4.8.2 洛杰斯蒂克(Logistic)模型要对马尔萨斯模型进行修

25、改，应进一步考虑哪些因素呢？人们发现在人口比较稀少，资源较丰富的条件下，人口增长较快，可以在短期内维持常数增长率。但当人口数量发展到一定水平后，会产生许多问题,如事物短缺，交通拥挤等，这又将导致人口增长率的减少，这种现象在某些动物种群的实验中也观察到。故可以假设，这反映了人口增长率随人口数量的增加而减少的现象，此时称为内禀增长率，为环境容纳量。按照这个修改的假设，就得到人口增长的洛杰斯蒂克(Logistic)模型 (4.22)对方程(4.22)求解得到 (4.23)从解的表达式(4.23)可以得到如下两点结论：(1) 其人口学含义是，不管开始时人口处于什么状态,当时间无限制增加时,人口总数都会趋于其环境容纳量。(2) 当时，,当时，其人口学含义是：当人口数量超过环境容纳量时，人口的数量将减少；而当人口数量小于环境容纳量时，人口数量将增加。图4-2和图4-3分别描述了人口总数及人口的增长率与时间的关系.图4-2反映了不管或,总有图4-3表明人口增长率的最大值在处。用1982年时我国人口统计作为初值，以36亿作为我国人口的容纳量，就可以得到表4-1中的预测值。从表4-1中的数据可以看出，作为短期预测，马尔萨斯模型与洛杰斯蒂克模型不相上下，但作为长期预报，洛杰斯蒂克模型比马尔萨斯模型要合理一些，