微分方程模型導(dǎo)彈跟蹤

導(dǎo)彈跟蹤問題常微分方程模型

問題某軍隊一導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處海面上有敵艇一艘以90km/h的速度向正東方向行駛。該基地立即發(fā)射導(dǎo)彈跟蹤追擊敵艇，導(dǎo)彈速度為450km/h。自動導(dǎo)航系統(tǒng)使導(dǎo)彈在任意時刻都能對準(zhǔn)敵艇。試問導(dǎo)彈在何時何處擊中敵艇？數(shù)學(xué)建模

微分方程建模的方法主要是依據(jù)守恒定律來建立等量關(guān)系式。對于這個問題，尋求等量關(guān)系是比較簡單的。設(shè)坐標(biāo)系如下圖所示，取導(dǎo)彈基地為原點(diǎn)(0,0)，x軸指向正東方，y軸指向正北方向。

P(x,y)

O

A(0,H)

x

y

M

當(dāng)t=0時，導(dǎo)彈位于O，敵艇位于點(diǎn)A(0,H)，其中H=120(km).設(shè)導(dǎo)彈在t時刻的位置為p(x(t),y(t)),由題意

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

P(x,y)

O

A(0,H)

x

M

y

另外在t時刻,敵艇位置應(yīng)為M(vet,H),其中ve=90(km/h)。由于導(dǎo)彈軌跡的切線方向必須指向敵艦，即直線PM的方向就是導(dǎo)彈軌跡上點(diǎn)p的切線方向，故有（3.2）或?qū)憺椋?.3）方程（3.1）、（3.3）連同初值條件x(0)=0,y(0)=0，構(gòu)成了一個關(guān)于時間變量t的一階微分方程組的初值問題.

（3.4）為了尋求x與y的關(guān)系，要設(shè)法消去變量t,由式（3.2）（3.2）得兩邊對t求導(dǎo)即有把式（3.1）（3.1）改寫為代入上式，就得到軌跡方程。這是一個二階非線性微分方程，加上初值條件，則得到導(dǎo)彈軌跡的數(shù)學(xué)模型(3.5)(3.6)(3.7)模型求解解法一：解析解法方程(3.5)可以降階,令記則（3.5）化為一階可分離變量方程即

兩邊積分可得：由初值條件（3.7）p|y=0=0得C=1,從而:上式通過分子有理化可改寫為兩式相加得到這樣我們又得到一個可分離的變量方程（3.8）兩端積分得到

利用x|y=0=0得到于是導(dǎo)彈軌跡方程為（3.9）設(shè)導(dǎo)彈擊中敵艇于B（L，H）P(x,y)

O

A(0,H)

x

y

MB(L,H)

*以y=H代入（3.9）式，得（3.10）而導(dǎo)彈擊中敵艇的時刻（3.11）將數(shù)據(jù)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分鐘后擊中敵艇。

能用解析方法求解非線性常微分方程固然不錯，但這樣的結(jié)果并不具有普遍性，因此最后還是需要用數(shù)值解方法進(jìn)行求解。下面我們再考慮這個問題的數(shù)值解，并與精確解作比較。

解法二：數(shù)值解法將初值問題（3.5）---（3.7）化為一階微分方程組

（3.12）（3.13）（3.14）取自變量y的步長為h=H/n，于是得分割點(diǎn)y0=0，y1=h，y2=2h，…，yn=nh=H下面介紹兩種近似算法來進(jìn)行數(shù)值處理。(1).Euler方法

Euler方法十分簡單，就是利用數(shù)值積分給出計算公式。對于第一個方程在分割區(qū)間[yk,yk+1]上進(jìn)行積分計算，有

第二個方程可同樣處理。設(shè)導(dǎo)彈到達(dá)(xk,yk)處的時刻為tk

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

O

A(0,H)

x

y

MB(L,H)

*通過迭代實(shí)現(xiàn)上面的算法，其中的相關(guān)數(shù)據(jù)為

H=120(km),ve=90(km/h),vw=450(km/h),λ=ve/vw設(shè)計程序時需要用到兩個數(shù)組x=(x(1),x(2),…,x(n+1))和p=(p(1),p(2),…,p(n+1))當(dāng)y=H，也就是y(n+1)=H時，導(dǎo)彈擊中艦艇。這時敵艦艇向東跑的距離約為

L=x(n+1)，所用的時間為T=L/ve。實(shí)現(xiàn)該算法的程序如下：daodangenzong1.m%daodangenzong1.mH=120;ve=90;vw=450;lamda=ve/vw;n=4;%將y的變化區(qū)間[0,H]進(jìn)行等分，可取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)%導(dǎo)彈擊中敵艦艇時，艦艇向東走的距離T=L/ve%導(dǎo)彈擊中敵艦艇時所用的時間右表是取n=4時的計算結(jié)果kykxkpk000013000.052601.50.123905.00.22412011.50.42此時：L≈11.5(km),T≈0.128(h)精確解是：L=25(km),T≈0.2778(h)結(jié)果不理想！再將n取得大一些計算。下表是對于不同的n值所對應(yīng)的計算結(jié)果。顯然，n越大（即h越?。?，結(jié)果就越精確。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)很接近了。但我們還可以用進(jìn)度更高的改進(jìn)Euler法進(jìn)行求解。(2).改進(jìn)的Euler方法（預(yù)估—校正法）

以一維情況為例，對問題Euler迭代格式是xk+1=xk+hf(xk,tk),其中h=△t,tk=t0+kh。Euler方法用的是左矩形求積公式，下面我們用梯形求積公式來改進(jìn)Eule迭代格式。

也就是其中的xk+1是未知的，我們做以下改動，設(shè)

且令于是對如下問題,可以寫出相應(yīng)的改進(jìn)Euler迭代格式

（3.12）（3.13）（3.14）編寫程序?qū)崿F(xiàn)以上改進(jìn)Euler法：daodangenzong2.m%daodangenzong2.mH=120;ve=90;vw=450;lamda=ve/vw;n=4;%將y的變化區(qū)間[0,H]進(jìn)行等分，可取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)%導(dǎo)彈擊中敵艦艇時，艦艇向東走的距離T=L/ve%導(dǎo)彈擊中敵艦艇時所用的時間將y的變化區(qū)間[0,H]4等分，計算結(jié)果如下：kykxkpk00001300.75000.05842603.50300.14223909.28270.2956412021.2781Inf與精確解L=25(km),T≈0.2778(h)相比，改進(jìn)Euler法收斂的更快。同樣在n=4時，Euler法的結(jié)果為

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

下表給出了不同等分下改進(jìn)Euler法的計算結(jié)果：

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法的計算結(jié)果：

與精確解L=25(km),T≈0.2778(h)相比，同樣等分下，改進(jìn)Euler法收斂的更快。在n=4時,下圖畫出了導(dǎo)彈軌跡由解析式所給出的精確曲線以及由Euler法和改進(jìn)的Euler法進(jìn)行數(shù)值計算的近似曲線。Euler法改進(jìn)Euler解析法可見在n=4時,數(shù)值解法還達(dá)不到精度要求,即不能擊中敵艇。下圖是在n=240時的計算結(jié)果,可見計算值與理論值相符的很好，即按照計算值發(fā)射導(dǎo)彈可以確保擊中敵艦艇。Euler法改進(jìn)Euler解析法解法三：仿真方法

如果建立微分方程很困難，或者微分方程很復(fù)雜而較難做出數(shù)值處理，常?？梢杂梅抡娣椒?

所謂仿真方法，顧名思義，指的是模仿真實(shí)時間行為和過程的方法。在這個具體問題中，就是一步步地模擬導(dǎo)彈追蹤敵挺的實(shí)際過程。而計算機(jī)仿真，則是在計算機(jī)上通過相應(yīng)的程序和軟件來實(shí)現(xiàn)對事件運(yùn)行的實(shí)際過程的模擬。設(shè)導(dǎo)彈和敵艇在初始時刻（t=0）分別位于P0(0,0)

和M0(0,H)，此時導(dǎo)彈指向M0。而在t=τ時，導(dǎo)彈的位置為P1(x1,y1

)，其中x1=0,y1=vwτ，敵艇的位置則為M1(veτ

,H).這時導(dǎo)彈沿P1M1方向飛行，P1M1的斜角為P0x

yM0P1M1θ1在t=2τ時，導(dǎo)彈的位置為P2(x2,y2

)，其中P2M2θ2M3M4此時敵艇位置為M2(2veτ

,H)導(dǎo)彈沿P2M2方向飛行.P0x

yM0P1M1θ1P2M2θ2M3M4以此方式，一般地，設(shè)t=kτ時，導(dǎo)彈位置為Pk(xk,yk

)敵艇的位置則為Mk(kveτ

,H)導(dǎo)彈將沿PkMk方向飛行,

那么，PkMk的斜角為

從而t=(k+1)τ時，導(dǎo)彈位置為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時，仿真停止；或者事先給定誤差界ε，當(dāng)yk+1-H<ε時,仿真停止.這時對于τ=0.1，0.05，0.005，0.001和0.0001，用仿真迭代格式（3.15）-（3.17）進(jìn)行計算，結(jié)果如下。τ0.10.010.0050.0010.0001L22.6822.2725.6725.0525.00T0.25190.28070.28520.27830.2783改進(jìn)Euler算法的計算結(jié)果n4812244896120240L21.2822.9723.564.2024.5524.7524.7924.88T0.2360.2550.2620.2690.2730.275