【数学建模】导弹发射追击问题的数学模型.doc

数学建模竞赛承 诺 书我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B中选择一项填写）： A 我们的队号为： 27 参赛队员：1. 唐路明 2. 季凯 3. 闻莺 指导教师或指导教师组负责人： 数模组 日期： 2009 年 8 月 11 日评阅编号（由评阅老师评阅前进行编号）：数学建模竞赛编 号 专 用 页评阅编号：评阅记录：评阅人评分备注导弹发射追击问题的数学模型摘要本文对导弹发射追击敌机问题进行了求解和计算机模拟，以微分方程为理论基础，根据题目要求，提出基本假设，建立合理的模型，并通过分析在给定不同速度条件下的轨迹方程，得到发射空对空和地对空两种导弹击毁敌机的条件。问题（1），建立微分方程模型，化二阶方程为一阶方程，从而得到导弹轨迹的解析表达式，发射该种空对空导弹击中敌机的的条件范围是 (0, )，为飞机速率与导弹速率之比。同时利用MATLAB70仿真，对导弹追踪敌机的过程进行了计算机检验和模拟，所得结果与所求相符。问题（2），首先，建立三维空间直角坐标系，在任意时刻t确定了导弹和飞机的空间位置坐标后，将导弹速度分解，再根据高度与水平距离比值不变的关系，将问题转化为二维平面直角坐标系上的追击问题。然后与问题（1）的处理相似，用差微分方法即可得导弹的轨迹，最后再对两者速度比值进行讨论后，得发射该种地对空导弹击毁敌机的的条件范围是(0,)。并最终再用MATLAB70进行检验分析，以证明模型的合理性。 问题（3），在一定范围内，对m,n进行合理取值，利用MATLAB70进行仿真模拟，确定导弹击毁飞机的最小速度。然后，根据实际情况，假设飞机能够在一定区域范围内发现导弹，忽略其加速的时间过程，对模型做进一步的改进。另一方面，对模型做更一般化的改进，对平面内任意两点的位置（即导弹和飞机相对位置的不固定性）进行分析讨论，得出导弹击毁飞机的轨迹方程，并对问题一的模型进行了模拟检验。最后，对模型的优缺点进行了分析，并指出其广泛的实际应用范围，以说明模型的现实重要意义。关键词 微分方程 MATLAB70仿真模拟 迭代法 椭圆区域一、问题的提出某边防导弹基地的雷达发现位于其正东 N 公里处有一家来犯敌机正欲逃往正北方向 M 公里处的安全区。该基地的I型空对空追踪导弹和II型地对空追踪导弹均可针对目标随时自动调节追踪方向，截击敌机。但敌机一旦进入安全区后，由于电子干扰作用，I型、II型导弹均将失去追踪目标，无法将敌机击毁。1 如此时恰有一架携有I型空对空追踪导弹的、与敌机处于同一飞行高度的巡航飞机在空中，基地即下令巡航飞机发射I型追踪导弹击毁敌机。试确定导弹追踪敌机的轨迹，并在适当的假定下给出发射该种导弹击毁敌机的条件；2 如此时在基地即发射II型地对空追踪导弹去击毁敌机，试确定导弹追踪敌机的轨迹，并在适当的假定下，及发射该种导弹击毁敌机的条件；3 若导弹的速度可在发射前根据需要设定，对于不同的N、M取值，编写计算机程序（语言不限），利用计算得到的数据说明怎样的发射速度可确保击毁敌机。二、问题的分析导弹追踪敌机，如要击毁成功，必须以足够大的速度在敌机未到达安全区时追上。本题即旨在通过建立模型，能够对导弹的跟踪轨迹方程和击毁敌机条件做出合理解答。再通过计算机，对其追踪过程进行合理的仿真模拟。对于问题（1），建立平面直角坐标系后，导弹在x,y方向的速度有=，=，联立方程组后，将二阶微分方程化为一阶微分方程，即可求解。接着对所求解进行MATLAB编程检验，并对追踪过程进行计算机检验拟合。对于问题（2），地对空的导弹在追击敌机的过程中，存在一个高度差的问题。因此必须建立三维空间直角坐标系来标示确定两者的位置。导弹在整个追击过程中，垂直于y轴的速度方向始终指向敌机，故始终有dz/dx=(dz/dt)/(dx/dt)=h/n等式成立，从而可将问题转化为二维平面内的追击问题。接着即是用与问题（1）中相似的微分方程的处理方法，求得y后，将其与m比较，可得结果。对于问题（3），应用MATLAB70进行编程求解，运用迭代法原理，对m,n等距取值，编程求解不同情况下导弹击毁敌机的最小发射速度。 三、基本假设1导弹与敌机的运动为质点运动。2导弹与敌机在发射和飞行过程中始终保持匀速运动（即无阻力）。3在整个追击过程中，区域内无其他物体的干扰。4敌机在飞行过程中，对导弹的追踪始终未察觉。5导弹与敌机的质量不计。6查阅资料表明，飞机飞行速度可取1000公里/每小时。7当基地监测到敌机时，立刻发射导弹。8在敌机未到达安全区前，导弹只要追上敌机，即可将敌机击毁。四、符号说明符号符号说明导弹速率飞机速率飞机速率与导弹速率的比飞机被监测到的点飞机与基地的水平距离M敌机安全的边界点飞机向正北方逃跑的安全距离飞机飞行高度飞机发现导弹的区域半径导弹与飞机的距离任意时刻t飞机在轴上的距离任意时刻t导弹在轴上的距离任意时刻t导弹在上的距离 飞机加速后的速度注：其余符号在文中出现时再予以说明。五、模型的建立与求解51 问题（1）的模型建立与求解 模型建立 建立直角坐标系如图一，设在t=0时刻导弹发现飞机，此时导弹的位置在，飞机的位置在,飞机以速度平行于y轴正向行驶，导弹以速度按指向飞机的方向发射飞行（>）.在任意时刻t导弹位于点，而飞机到达Q（n, t）点，直线PQ与导弹轨迹（图一中曲线）相切，切线与x轴正向夹角为. 导弹在x,y方向的速度分别为，由直角三角形PQR写出sin（）和cos（）的表达式，得到微分方程 （式1）初始条件为 ， 图一模型求解 要想得到导弹拦截飞机的精确时间和位置，必须对问题一的求解作进一步的分析，设法得到某种形式的解析解。 将方程(1)的两式相除，消去得到 （式2）对x求导得 （式3）为消去式中的dt，利用曲线弧长s的微分，导弹的速度，以及微分关系，得 （式4）这是关于y(x)的二阶微分方程，若令，可化为的一阶微分方程， （式5）依题意其初始条件为. （式6）方程（式5）为可分离变量方程，容易求得在（式6）下的解为 （式7）或 （式8）由（式7）和（式8）得 （式9）对（式9）积分并注意到，得 （式10）此即为导弹轨迹的解析表达式，并用MATLAB编程进行检验，结果与所求相同。程序见附录一当=n时， 有， （式11）1）当k>1（即飞机的速率大于导弹的速率）时，方程的解为显然导弹不可能击中飞机,符合实际情况;2）当k=1（即飞机的速率等于导弹的速率）时，方程的解为 显然,表明此时导弹不可能击中飞机,符合实际情况;3）当k<1（即飞机的速率小于导弹的速率）时，方程的解为此时,当<m时,即 <m得 （式12）如图二所示 图二又 >0 故 的取值范围是(0, ). 模型检验 利用MATLAB编程仿真检验导弹轨迹与所求得解析解y的拟合程度，结果见图三。程序见附录二。 注：（1）经查阅资料，飞机的飞行速度约为1000公里/每小时； （2）在给定速度1000公里/每小时的情况下，取一较小值t的时间间隔为t=0.000005小时，则所求步长为d=0.000005*1000=0.005公里。运行结果显示，当d<=0.005时，导弹运动轨迹越过飞机航线，继续追踪则出现折转。且每隔0.000005秒，导弹轨迹都折转一次，以飞机航线为轴，成“Z”字型路线前进。因此，在程序中将t的值累计叠加，直至当d>0.005时，导弹运行轨迹才恢复正常，不发生折转。而为符合实际情况，不使d过大，在d<=0.01时，即认为导弹击中飞机。 图三 52 问题（2）的模型建立与求解 模型建立 建立空间三维直角坐标系如图三。设在t=0时刻，地对空的导弹发现了敌机，此时导弹在原点的位置，而飞机的位置在N（n,0,h）点（飞机高度为h），飞机以速度沿-y方向即XOY平面内正北方向行驶，导弹以速度按指向飞机的方向发射飞行。在任意时刻t，飞机位于R（n, ,h）点（），导弹位于P(,)点。此时沿垂直于轴作切面,得图中阴影部分。导弹P在阴影面上的投影点的坐标为P(,)。将t时刻导弹速度沿平行于y轴和垂直于y轴两个方向分解,其中平行于y轴的速度为（为速度方向与阴影面的夹角），方向为。由t的任意性可知，导弹速度在垂直于y轴方向上的速度方向始终不变，即dz/dx=(dz/dt)/(dx/dt)=h/n。 图四 由于导弹在XOZ平面内的速度方向始终不变，则导弹在XOZ平面内的运动轨迹也始终不变，即对于任意导弹位置P（x,y,z）恒有x/z=n/h。又x/z=n/h所表示的几何意义是一个如图五中所示斜平面，且导弹在该平面内运动。此时，导弹在三维空间内的运动轨迹即可转化为问题一中的二维平面追击问题，如图六所示。 图五 图六模型求解 由图六可得方程组 （式13）由问题一中结果得 （式14） 代入 ，得 （式15） 又 ，得 （式16）则导弹的轨迹方程为 （式17）当时， （式18）1）当k>1（即飞机的速率大于导弹的速率）时，方程的解为显然导弹不可能击中飞机,符合实际情况;2）当k=1（即飞机的速率等于导弹的速率）时，方程的解为 显然,表明此时导弹不可能击中飞机,符合实际情况;3）当k<1（即飞机的速率小于导弹的速率）时，方程的解为此时,当>m时,即 >m得 （式19） （如图七所示） 图七又 故 的取值范围是(0, ). 53 问题3的求解 根据题目信息所得，当给定N，M，值时，计算出导弹可击毁飞机的最小速度。因此，可应用MATLAB70进行过程仿真。算法如下：一、 确定飞机与基地的位置，设为，二、 输入m,n,三、 利用循环选取大量的，对每一个，给定时间间隔为，计算每一点在各个时刻的坐标。设某点在t时刻的坐标为，其中 四、 取足够小的时结束算法。五、 对每一个点，连接它 各个时刻的位置，即得所求运动的轨迹。注：考虑到计算机运行的效率，在选取时，可以分为两个步骤。先可以确定一个小范围（如在1的范围内），然后再将小范围细分，这样就可以减少大量的运算时间，并且得到的数据更加精确。所得结果如表一所示。（程序见附录三）NM20406080100201.6167e+0031.2788e+0031.1791e+0031.1317e+0031.1041e+003402.4178e+0031.6165e+0031.3862e+0031.2798e+0031.2190e+003603.3118e+0031.9990e+0031.6170e+0031.4421e+0031.3433e+003804.2504e+0032.4139e+0031.8677e+0031.6172e+0031.4763e+0031005.2123e+0032.8511e+0032.1345e+0031.8036e+0031.6174e+003表一六、模型的改进与优化改进一在上述对整个题目的求解过程中，假设飞机对导弹始终未察觉。而在实际情况中，飞机在逃亡时，也会勘察周边环境，根据具体情况对飞行路线和速度大小进行改变。故对第一问中所建模型的改进，同样也适用于本题的第二问。在此，仅对于问题一的模型改进进行阐述。模型建立 假设飞机在逃亡的过程中，若勘察到在它周围半径为r的范围内有追踪导弹，则会在最短时间内立即调整，调整后速度大小为，方向仍然不变。进而设导弹在刚好进入飞机半径为r的范围时所经过的时间为t，导弹的坐标为P（,），飞机的坐标Q（n，t）；则 在时刻t之前，导弹轨迹方程是 (式20) 在时刻t之后，设 ，此时可建立以P点为原点，正东方向为x方向，正北方向为y方向的直角坐标系，如下图所示，则在飞机改变速度后，导弹的运动轨迹为图中实线所示。图八此时（在时刻t后），导弹的轨迹方程是 (式21)其中 整理得 (式22) 故最后所得结论为（1）飞机未发现导弹前，即时刻t之前导弹轨迹方程为 （2）飞机发现导弹后，即时刻t之后导弹轨迹方程为其中， 模型的计算机模拟 利用MATLAB对导弹轨迹方程进行模拟，结果如图九所示。（程序见附录四）图九改进二将此题二维模型更一般处理，即设平面内有任意两点，分别为飞机和导弹的位置，建立如下图所示的坐标系。在任意t时刻，导弹P的位置为（，），运动方向始终指向飞机Q，速度大小为；飞机Q的位置为（，），即运动方向始终与x轴夹角为。且和大小均为常数，。出于一般性考虑，认为可以变化，则t时刻飞机的位置为： (式23)而导弹在追击过程中满足：， (式24)由（24）式可解出导弹的运动微分方程： (式25)可以由解微分方程组求出，从而得到导弹的运动轨迹。在此建立以飞机为原点，飞机运动方向为y轴正方向的直角坐标系，即相当于取，则由（23）式和（25）式消去时间t得到导弹的运动轨迹方程 （a，b为与初始量有关的参数）由于计算数值变化量过大，故不能再进一步自行讨论；根据参考文献【4】中P217221中计算模拟结论得到：该平面存在一个范围内S，若导弹的初始位置在S内，则导弹可以在飞机进入安全位置前击中它；若导弹的初始位置在S以外，则飞机可以安全返回安全位置，导弹就再也不能击中飞机。如下图所示：该图S有以下几个性质：（1）S图关于y轴对称，也有可能是一个椭圆；（2）；（3）因为导弹在图形S边界上任意一点出发，都是在点M击中飞机，即飞机所走的路程相同，则说明追击过程所花的时间相同。将此模型对问题一中的模型进行检验。因为此图形和椭圆的相似度很高，故把次边际认为是椭圆。则椭圆的中心为（0，m），焦点为（0，0）和（0，2m）。得到椭圆的方程为：取y=0时，有： ，又0<k<1化简，可得：由问题一中的导弹坐标为（-n，0），则：即：，与一式相同，则证明问题一中模型的合理性。七、模型的评价与推广 本文的模型主要建立在差微分方程求解的基础上，简单易懂，操作性强。在对击毁条件的求解过程中，模型对导弹与飞机的速率进行了详细讨论。并根据飞机是否察觉导弹的基本前提，对模型进行了切实有效的改进。在每次求解过后，都利用MATLAB70进行编程检验，确保结果更加真实可靠。所编程序运行效率较高，最长亦可在两分钟内得出准确结果。 但是，在实际情况中，不可能不考虑阻力、重力等因素；飞机在改变速度时，也必然有一个时间上的加速过程。 此类追击问题的模型建立，具有广泛的应用范围和深远的现实意义。在现实海陆空领域防卫中，敌我双方均可利用该模型来进行防御与攻击战略的调整。此外，该模型亦可有效适用于走私、海上缉毒、航母与护航舰搜寻人质问题，潜艇决策等侦察活动中，在生产与军事上有较好的实际意义。八、参考文献【1】 姜启源，邢文训等大学数学实验北京：清华大学出版社，2005【2】 岂兴明，王占富等MATLAB 7.0程序设计快速入门北京：人民邮电出版社，2009【3】 白其峥主编数学建模案例分析北京：海洋出版社，2000【4】 李继成主编数学实验北京：高等教育出版社，2006【5】 邹秀芬，陈绍林等数值计算方法学习指导书武汉：武汉大学出版社，2008【6】 韩中庚主编数学建模竞赛获奖论文精选与点评北京：科学出版社，2007【7】 熊启才主编数学模型方法及应用重庆：重庆大学出版社，20053【8】 王树禾编著数学模型基础合肥：中国科学技术大学出版社，1996附录：附录一 对于问题一中解析解y的检验程序 >> a=dsolve('Dy=0.5*(nk)/(n-x)k-(n-x)k/(nk)','y(0)=0','x') a = -1/2*nk*(n-x)(-k+1)/(-k+1)+1/2*n(-k)*(n-x)(k+1)/(k+1)-n*k/(k2-1) >> pretty(a) k (-k + 1) (-k) (k + 1) n (n - x) n (n - x) n k - 1/2 - + 1/2 - - - -k + 1 k + 1 2 k - 1>> 附录二 对于问题一中导弹追击飞机的检验程序v=1000;%定义速度为1000公里每小时dt=0.000005;%时间间隔为0.000005秒x=0 100;%定义初始的X坐标位置y=0 0;%定义初始的Y坐标位置t=0;M=70;%飞机安全的临界值for i=1:2 plot(x(i),y(i),'.'),hold onendd=100;%飞机与导弹的初始距离while(d>0.01)%当飞机与导弹的距离小于等于0.01时，结束循环 t=t+dt;%时间的叠加 d=sqrt(x(2)-x(1)2+(y(2)-y(1)2);%飞机与导弹的距离 x(1)=x(1)+2*v*dt*(x(2)-x(1)/d;%导弹X的位置变化 y(1)=y(1)+2*v*dt*(y(2)-y(1)/d;%导弹Y的位置变化 y(2)=y(2)+v*dt;%飞机Y的位置变化 x(2)=x(2);%飞机X位置没有变化 for i=1:2 plot(x(i),y(i),'.'),hold on endendif y(2)>M %判断飞机是否安全 fprintf('导弹无法击毁飞机')else fprintf('导弹可以击毁飞机')endtxlabel('x')ylabel('y')title('导弹追击飞机的轨迹的检验')hold onx=0:1:100;y=(100/3).*(100-x)./100).(3/2)-100.*(100-x)./100).(1/2)+200/3;plot(x,y,'r+');保存为dd.m在Command Windou中输入 dd 显示如下：>> dd导弹可以击毁飞机t = 0.0667附录三 问题三的求解程序input('该程序可以解决导弹追踪问题中导弹最小速度的确定，M是安全的边界值，N是基地发现敌机时之间的距离')input('飞机的速度我们取了1000公里每小时，您也可以自行调整并且我们还可以为您画出具体的追踪轨迹，请按提示进行操作')M=input('输入M：')N=input('输入N：')v=1000;%定义速度为1000公里每小时dt=0.00005;%时间间隔为0.000005小时x=0 N;%定义初始的X坐标位置y=0 0;%定义初始的Y坐标位置d=N;%飞机与导弹的初始距离for k=1001:1:10000 while(d>0.1)%当飞机与导弹的距离小于等于0.1时，结束循环 d=sqrt(x(2)-x(1)2+(y(2)-y(1)2);%飞机与导弹的距离 x(1)=x(1)+k*dt*(x(2)-x(1)/d;%导弹X的位置变化 y(1)=y(1)+k*dt*(y(2)-y(1)/d;%导弹Y的位置变化 y(2)=y(2)+v*dt;%飞机Y的位置变化 x(2)=x(2);%飞机X位置没有变化 end if y(2)<M break; end x(1)=0; y(1)=0; y(2)=0; d=N;endfor p=k-1:0.0001:kv=1000;%定义速度为1000公里每小时dt=0.00005;%时间间隔为0.000005小时x=0 N;%定义初始的X坐标位置y=0 0;%定义初始的Y坐标位置d=N;%飞机与导弹的初始距离 while(d>0.1)%当飞机与导弹的距离小于等于0.01时，结束循环 d=sqrt(x(2)-x(1)2+(y(2)-y(1)2);%飞机与导弹的距离 x(1)=x(1)+p*dt*(x(2)-x(1)/d;%导弹X的位置变化 y(1)=y(1)+p*dt*(y(2)-y(1)/d;%导弹Y的位置变化 y(2)=y(2)+v*dt;%飞机Y的位置变化 x(2)=x(2);%飞机X位置没有变化 end if y(2)<M break; end x(1)=0; y(1)=0; y(2)=0; d=N;end disp('当发射速度大于')pdisp('可确保导弹击中飞机')x=0 N;%定义初试X的坐标位置y=0 0;%定义初试Y的坐标位置t=0;for i=1:2 plot(x(i),y(i),'.'),hold onendd=N;%飞机与导弹的初始距离while(d>0.1)%当飞机与导弹的距离小于等于0.01公里时，循环结束 t=t+dt; d=sqrt(x(2)-x(1)2+(y(2)-y(1)2); x(1)=x(1)+p*dt*(x(2)-x(1)/d; y(1)=y(1)+p*dt*(y(2)-y(1)/d; y(2)=y(2)+v*dt; x(2)=x(2); for i=1:2 plot(x(i),y(i),'.'),hold on endendxlabel('x')ylabel('y')title('导弹追踪飞机的轨迹')附录四 问题一的模型改进计算机检验程序v=1000;%定义速度为1000公里每小时dt=0.000005;%时间间隔为0.000005秒x=0 100;%定义初试X的坐标位置y=0 0;%定义初试Y的坐标位置t=0;M=70;%定义飞机安全的临界值for i=1:2 plot(x(i),y(i),'.'),hold onendd=100;%飞机与导弹的初始距离while(d>0.01)%当飞机与导弹的距离小于等于0.01公里时，循环结束 t=t+dt; if(d>5)%当飞机与导弹的距离大于5公里时，飞机没有察觉有导弹追踪他 d=sqrt(x(2)-x(1)2+(y(2)-y(1)2); x(1)=x(1)+2*v*dt*(x(2)-x(1)/d; y(1)=y(1)+2*v*dt*(y(2)-y(1)/d; y(2)=y(2)+v*dt; x(2)=x(2); for i=1:2 plot(x(i),y(i),'.'),hold on end else %当飞机与导弹的距离小于等于5公里时，飞机加速 d=sqrt(x(2)-x(1)2+(y(2)-y(1)2); x(1)=x(1)+2*v*dt*(x(2)-x(1)/d; y(1)=y(1)+2*v*dt*(y(2)-y(1)/d; y(2)=y(2)+1.5*v*dt; x(2)=x(2); for i=1:2 plot(x(i),y(i),'.'),hold on end endendif y(2)>M fprintf('导弹无法击毁飞机')else fprintf('导弹可以击毁飞机')endtxlabel('x')ylabel('y')title('导弹追踪飞机的轨迹')保存为aa.m在Command Window 里输入aa 显示如下：cc导弹无法击毁飞机t = 0.0716>>