導彈跟蹤問題常微分方程模型

問題某軍隊一導彈基地發現正北方向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