微分方程模型导弹跟踪

导弹跟踪问题常微分方程模型

问题某军队一导弹基地发现正北方向120km处海面上有敌艇一艘以90km/h的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹速度为450km/h。自动导航系统使导弹在任意时刻都能对准敌艇。试问导弹在何时何处击中敌艇？数学建模

微分方程建模的方法主要是依据守恒定律来建立等量关系式。对于这个问题，寻求等量关系是比较简单的。设坐标系如下图所示，取导弹基地为原点(0,0)，x轴指向正东方，y轴指向正北方向。

P(x,y)

O

A(0,H)

x

y

M

当t=0时，导弹位于O，敌艇位于点A(0,H)，其中H=120(km).设导弹在t时刻的位置为p(x(t),y(t)),由题意

（3.1）其中vw=450(km/h)。y

P(x,y)

O

A(0,H)

x

M

y

另外在t时刻,敌艇位置应为M(vet,H),其中ve=90(km/h)。由于导弹轨迹的切线方向必须指向敌舰，即直线PM的方向就是导弹轨迹上点p的切线方向，故有（3.2）或写为（3.3）方程（3.1）、（3.3）连同初值条件x(0)=0,y(0)=0，构成了一个关于时间变量t的一阶微分方程组的初值问题.

（3.4）为了寻求x与y的关系，要设法消去变量t,由式（3.2）（3.2）得两边对t求导即有把式（3.1）（3.1）改写为代入上式，就得到轨迹方程。这是一个二阶非线性微分方程，加上初值条件，则得到导弹轨迹的数学模型(3.5)(3.6)(3.7)模型求解解法一：解析解法方程(3.5)可以降阶,令记则（3.5）化为一阶可分离变量方程即

两边积分可得：由初值条件（3.7）p|y=0=0得C=1,从而:上式通过分子有理化可改写为两式相加得到这样我们又得到一个可分离的变量方程（3.8）两端积分得到

利用x|y=0=0得到于是导弹轨迹方程为（3.9）设导弹击中敌艇于B（L，H）P(x,y)

O

A(0,H)

x

y

MB(L,H)

*以y=H代入（3.9）式，得（3.10）而导弹击中敌艇的时刻（3.11）将数据H=120(km),ve=90(km/h),vw=450(km/h)代入(3.10)、(3.11)式，得到L=25(km),T≈0.2778(h)=13分钟即在东面25公里处，13分钟后击中敌艇。

能用解析方法求解非线性常微分方程固然不错，但这样的结果并不具有普遍性，因此最后还是需要用数值解方法进行求解。下面我们再考虑这个问题的数值解，并与精确解作比较。

解法二：数值解法将初值问题（3.5）---（3.7）化为一阶微分方程组

（3.12）（3.13）（3.14）取自变量y的步长为h=H/n，于是得分割点y0=0，y1=h，y2=2h，…，yn=nh=H下面介绍两种近似算法来进行数值处理。(1).Euler方法

Euler方法十分简单，就是利用数值积分给出计算公式。对于第一个方程在分割区间[yk,yk+1]上进行积分计算，有

第二个方程可同样处理。设导弹到达(xk,yk)处的时刻为tk

,那么得到计算的迭代格式。3.153.173.16P(x,y)

O

A(0,H)

x

y

MB(L,H)

*通过迭代实现上面的算法，其中的相关数据为

H=120(km),ve=90(km/h),vw=450(km/h),λ=ve/vw设计程序时需要用到两个数组x=(x(1),x(2),…,x(n+1))和p=(p(1),p(2),…,p(n+1))当y=H，也就是y(n+1)=H时，导弹击中舰艇。这时敌舰艇向东跑的距离约为

L=x(n+1)，所用的时间为T=L/ve。实现该算法的程序如下：daodangenzong1.m%daodangenzong1.mH=120;ve=90;vw=450;lamda=ve/vw;n=4;%将y的变化区间[0,H]进行等分，可取n=4,8,…,240h=H/n;x(1)=0;p(1)=0;y=0:h:H;fork=1:nx(k+1)=x(k)+h*p(k);p(k+1)=p(k)+h*lamda*sqrt(1+p(k)^2)/(H-y(k));endL=x(n+1)%导弹击中敌舰艇时，舰艇向东走的距离T=L/ve%导弹击中敌舰艇时所用的时间右表是取n=4时的计算结果kykxkpk000013000.052601.50.123905.00.22412011.50.42此时：L≈11.5(km),T≈0.128(h)精确解是：L=25(km),T≈0.2778(h)结果不理想！再将n取得大一些计算。下表是对于不同的n值所对应的计算结果。显然，n越大（即h越小），结果就越精确。n4812244896120240L11.5215.9617.9720.2522.2523.3323.5824.15T0.1280.1770.2000.2880.2470.2590.2620.268此时的近似解：L≈24.15(km),T≈0.268(h)与精确解L=25(km),T≈0.2778(h)很接近了。但我们还可以用进度更高的改进Euler法进行求解。(2).改进的Euler方法（预估—校正法）

以一维情况为例，对问题Euler迭代格式是xk+1=xk+hf(xk,tk),其中h=△t,tk=t0+kh。Euler方法用的是左矩形求积公式，下面我们用梯形求积公式来改进Eule迭代格式。

也就是其中的xk+1是未知的，我们做以下改动，设

且令于是对如下问题,可以写出相应的改进Euler迭代格式

（3.12）（3.13）（3.14）编写程序实现以上改进Euler法：daodangenzong2.m%daodangenzong2.mH=120;ve=90;vw=450;lamda=ve/vw;n=4;%将y的变化区间[0,H]进行等分，可取n=4,8,…,240h=H/n;x(1)=0;p(1)=0;y=0:h:H;fork=1:nx1(k+1)=x(k)+h*p(k);p1(k+1)=p(k)+h*lamda*sqrt(1+p(k)^2)/(H-y(k));x(k+1)=1/2*(x1(k+1)+x(k)+h*p1(k+1))p(k+1)=1/2*(p1(k+1)+p(k)+h*lamda*sqrt(1+p1(k+1)^2)/(H-y(k+1)))endL=x(n+1)%导弹击中敌舰艇时，舰艇向东走的距离T=L/ve%导弹击中敌舰艇时所用的时间将y的变化区间[0,H]4等分，计算结果如下：kykxkpk00001300.75000.05842603.50300.14223909.28270.2956412021.2781Inf与精确解L=25(km),T≈0.2778(h)相比，改进Euler法收敛的更快。同样在n=4时，Euler法的结果为

L≈11.5(km),T≈0.128(h)此时L≈x4=21.2781,T≈L/ve=0.2364

下表给出了不同等分下改进Euler法的计算结果：

n4812244896120240L21.2822.9723.564.2024.5524.7524.7924.88T0.2360.2550.2620.2690.2730.2750.2760.277n4812244896120240L11.5215.9617.9720.2522.2523.3323.5824.15T0.1280.1770.2000.2880.2470.2590.2620.268下表是Euler法的计算结果：

与精确解L=25(km),T≈0.2778(h)相比，同样等分下，改进Euler法收敛的更快。在n=4时,下图画出了导弹轨迹由解析式所给出的精确曲线以及由Euler法和改进的Euler法进行数值计算的近似曲线。Euler法改进Euler解析法可见在n=4时,数值解法还达不到精度要求,即不能击中敌艇。下图是在n=240时的计算结果,可见计算值与理论值相符的很好，即按照计算值发射导弹可以确保击中敌舰艇。Euler法改进Euler解析法解法三：仿真方法

如果建立微分方程很困难，或者微分方程很复杂而较难做出数值处理，常常可以用仿真方法.

所谓仿真方法，顾名思义，指的是模仿真实时间行为和过程的方法。在这个具体问题中，就是一步步地模拟导弹追踪敌挺的实际过程。而计算机仿真，则是在计算机上通过相应的程序和软件来实现对事件运行的实际过程的模拟。设导弹和敌艇在初始时刻（t=0）分别位于P0(0,0)

和M0(0,H)，此时导弹指向M0。而在t=τ时，导弹的位置为P1(x1,y1

)，其中x1=0,y1=vwτ，敌艇的位置则为M1(veτ

,H).这时导弹沿P1M1方向飞行，P1M1的斜角为P0x

yM0P1M1θ1在t=2τ时，导弹的位置为P2(x2,y2

)，其中P2M2θ2M3M4此时敌艇位置为M2(2veτ

,H)导弹沿P2M2方向飞行.P0x

yM0P1M1θ1P2M2θ2M3M4以此方式，一般地，设t=kτ时，导弹位置为Pk(xk,yk

)敌艇的位置则为Mk(kveτ

,H)导弹将沿PkMk方向飞行,

那么，PkMk的斜角为

从而t=(k+1)τ时，导弹位置为Pk+1(xk+1,yk+1

)，其中

（3.15）（3.16）（3.17）（3.18）而敌艇位置为

Mk+1((k+1)veτ

,H)。计算直至yk<H,yk+1≥H时，仿真停止；或者事先给定误差界ε，当yk+1-H<ε时,仿真停止.这时对于τ=0.1，0.05，0.005，0.001和0.0001，用仿真迭代格式（3.15）-（3.17）进行计算，结果如下。τ0.10.010.0050.0010.0001L22.6822.2725.6725.0525.00T0.25190.28070.28520.27830.2783改进Euler算法的计算结果n4812244896120240L21.2822.9723.564.2024.5524.7524.7924.88T0.2360.2550.2620.2690.2730.275