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

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

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

微分方程建模的方法主要是依據(jù)守恒定律來(lái)建立等量關(guān)系式。對(duì)于這個(gè)問(wèn)題，尋求等量關(guān)系是比較簡(jiǎ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時(shí)，導(dǎo)彈位于O，敵艇位于點(diǎn)A(0,H)，其中H=120(km).設(shè)導(dǎo)彈在t時(shí)刻的位置為p(x(t),y(t)),由題意

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

P(x,y)

O

A(0,H)

x

M

y

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

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

兩邊積分可得：由初值條件（3.7）p|y=0=0得C=1,從而:上式通過(guò)分子有理化可改寫(xiě)為兩式相加得到這樣我們又得到一個(gè)可分離的變量方程（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)彈擊中敵艇的時(shí)刻（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分鐘后擊中敵艇。

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

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

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

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

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

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

O

A(0,H)

x

y

MB(L,H)

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

H=120(km),ve=90(km/h),vw=450(km/h),λ=ve/vw設(shè)計(jì)程序時(shí)需要用到兩個(gè)數(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時(shí)，導(dǎo)彈擊中艦艇。這時(shí)敵艦艇向東跑的距離約為

L=x(n+1)，所用的時(shí)間為T(mén)=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)彈擊中敵艦艇時(shí)，艦艇向東走的距離T=L/ve%導(dǎo)彈擊中敵艦艇時(shí)所用的時(shí)間右表是取n=4時(shí)的計(jì)算結(jié)果kykxkpk000013000.052601.50.123905.00.22412011.50.42此時(shí)：L≈11.5(km),T≈0.128(h)精確解是：L=25(km),T≈0.2778(h)結(jié)果不理想！再將n取得大一些計(jì)算。下表是對(duì)于不同的n值所對(duì)應(yīng)的計(jì)算結(jié)果。顯然，n越大（即h越小），結(jié)果就越精確。n4812244896120240L11.5215.9617.9720.2522.2523.3323.5824.15T0.1280.1770.2000.2880.2470.2590.2620.268此時(shí)的近似解：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ù)估—校正法）

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

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

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

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

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

下表給出了不同等分下改進(jìn)Euler法的計(jì)算結(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法的計(jì)算結(jié)果：

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

如果建立微分方程很困難，或者微分方程很復(fù)雜而較難做出數(shù)值處理，常常可以用仿真方法.

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

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

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

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

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

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

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

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

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

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

那么，PkMk的斜角為

從而t=(k+1)τ時(shí)，導(dǎo)彈位置為Pk+1(xk+1,yk+1

)，其中

（3.15）（3.16）（3.17）（3.18）而敵艇位置為

Mk+1((k+1)veτ

,H)。計(jì)算直至yk<H,yk+1≥H時(shí)，仿真停止；或者事先給定誤差界ε，當(dāng)yk+1-H<ε時(shí),仿真停止.這時(shí)對(duì)于τ=0.1，0.05，0.005，0.001和0.0001，用仿真迭代格式（3.15）-（3.17）進(jìn)行計(jì)算，結(jié)果如下。τ0.10.010.0050.0010.0001L22.6822.2725.6725.0525.00T0.25190.28070.28520.27830.2783改進(jìn)Euler算法的計(jì)算結(jié)果n4812244896120240L21.2822.9723.564.2024.5524.7524.7924.88T0.2360.2550.2620.2690.2730.275