魚雷擊艦問題數(shù)模論文

姓名：班級：聯(lián)系電話：姓名：班級：聯(lián)系電話：姓名：班級：聯(lián)系電話：姓名班級負(fù)責(zé)項目學(xué)習(xí)課程（注：關(guān)于數(shù)學(xué)、物理、經(jīng)濟(jì)、編程的課程）班級排名自動化創(chuàng)新班建模10電子信息工程1班論文10自動化1班程序魚雷擊艦問題數(shù)模論文谷小龍；陳磊；鄧博文（機(jī)電學(xué)院自動化創(chuàng)新班，機(jī)電學(xué)院自動化班，機(jī)電學(xué)院電子班）摘要：論文討論的是魚雷追蹤艦艇的運動軌跡。首先，我們應(yīng)對魚雷的追蹤軌跡進(jìn)行分析，判斷屬于哪種類別函數(shù)。其次，建立微分方程模型，是一個二階方程,通過降階法化為一階方程，然后用分離變量法求解。得到所求的方程表達(dá)式，再結(jié)合魚雷速率和艦艇速率，因而得到魚雷追蹤艦艇的軌跡方程。然后，在假定忽略萬有引力與水的阻力的條件下，建立數(shù)學(xué)模型。利用matlab數(shù)學(xué)軟件編程求解參數(shù)方程，得到敵艦的航行距離約為0.667海里，進(jìn)一步求得時間約等于1.48min。針對本題我們建立的微分方程模型符合題目的要求，并且我們通過解析法和參數(shù)方程MATLAB求值的方法得到了基本相同的結(jié)果（結(jié)果的誤差率不超過1%），并在內(nèi)容中附上了相應(yīng)的曲線圖與原始數(shù)據(jù)。這同時也論證了我們的模型的正確性、合理性。針對本問題我們還對其進(jìn)行了拓展，假設(shè)敵艦改變了速度與方向，與魚雷方向成一夾角逃跑的情形，通過仿真的方法建立了合理的模型，使得該問題更貼近于現(xiàn)實的實際情況，有一定的應(yīng)用意義。最后，我們對模型進(jìn)行了評價與改進(jìn)，并對模型中的參數(shù)進(jìn)行了分析，最終得出了結(jié)論“魚雷在1.48min之后將擊中敵艦”。另外，在附錄中附有matlab的源程序與改進(jìn)的原始數(shù)據(jù)。關(guān)鍵詞：微分方程，魚雷速率，艦艇速率，matlab，參數(shù)方程數(shù)值法仿真拓展一、問題的重述魚雷擊艦問題(微分方程)如圖所示,一敵艦在某海域內(nèi)沿正北方向航行時,我方戰(zhàn)艦恰位于敵艦正西方1海里處,我艦向敵艦發(fā)射制導(dǎo)魚雷,敵艦速度為0.45海里/分鐘,魚雷速度為敵艦速度的2倍,魚雷的運行方向始終指向敵艦,試問敵艦航行多遠(yuǎn)時將被擊中?二、問題的分析根據(jù)題意，畫出示意圖如上。由題意可知，敵艦的速度方向與大小不變，恒為0.45海里每分鐘，魚雷的速度大小也恒定為敵艦的2倍，即0.9海里每分鐘。由魚雷的運行方向始終指向敵艦可知，魚雷的運行方向隨著敵艦的位置不同而發(fā)生著改變，這就是本題的難點所在。對于非線性曲線的問題，用常規(guī)的方法難以模擬出來，所以我們采用引入微積分，用微積分的表達(dá)式表示圖示中變量的方程。微積分可以使模擬過程簡化，而且易于表述出來。如上圖中取曲線上任意一點P（x（t），y（t）），敵艦位于Q（l，vt）由于魚雷頭始終對準(zhǔn)敵艦，故此時直線PQ就是魚雷的軌跡曲線弧OP在點P處的切線，即有y=，即=（1—x）y+y又根據(jù)題意，弧OP的長度的2倍，即dx=2由以上兩式消去t，整理得模型：（1—x）y=。這就是用微積分方程建立的模型，表述簡單易于理解。但是最后的模型是一個二階微分方程，求解起來有一定的難度，過程有點復(fù)雜。所以除了用解析法將方程直接求解出來，再將數(shù)據(jù)代入求解之外，我們還引用了MATLAB軟件直接求解，過程簡單、直接將結(jié)果表述出來，但結(jié)果會存在一定的誤差。三、基本假設(shè)假設(shè)：魚雷的質(zhì)量忽略不記，在魚雷及敵艦運動過程中忽略水流和氣流向?qū)λ鼈冞\動軌跡的影響。四、定義符號說明定義說明：微積分：微積分（Calculus）是高等數(shù)學(xué)中研究函數(shù)的微分、積分以及有關(guān)概念和應(yīng)用的數(shù)學(xué)分支。微積分學(xué)是微分學(xué)和積分學(xué)的總稱。它是一種數(shù)學(xué)思想，‘無限細(xì)分’就是微分，‘無限求和’就是積分。速度，時間，距離公式：當(dāng)時間一定距離與速度成正比。解析法：數(shù)學(xué)中用解析式表示函數(shù)的方法。參數(shù)方程：在給定的平面直角坐標(biāo)系中，如果曲線上任意一點的坐標(biāo)x，y都是某個變數(shù)t的函數(shù)x=f(t),y=φ(t)——(1)；且對于t的每一個允許值，由方程組(1)所確定的點m(x，y)都在這條曲線上，那么方程組(1)稱為這條曲線的參數(shù)方程，聯(lián)系x、y之間關(guān)系的變數(shù)稱為參變數(shù)，簡稱參數(shù)。類似地，也有曲線的極坐標(biāo)參數(shù)方程ρ=f(t),θ=g(t)。數(shù)值法：數(shù)值法主要是指有限元法,有限元法大多是在解析法FGM模型的基礎(chǔ)上,在不同尺度上進(jìn)行有限元離散,離散單元尺度不同,進(jìn)行有限元計算時要滿足的連續(xù)性條件不同,預(yù)測結(jié)果的精確度就不同符號說明：t-------任意時刻P（x（t），y（t））------t時刻魚雷位置Q（l，vt）--------t時刻敵艦位置PQ--------曲線上P點上的切線A(1,0)--------敵艦最初位置v----------敵艦的速度X----------t時刻魚雷的橫坐標(biāo)y----------t時刻魚雷的縱坐標(biāo)y----------y的一階求導(dǎo)y----------y的二階求導(dǎo)五、模型的建立圖示（建立坐標(biāo)系）假設(shè)t時刻魚雷的位置為P（x（t），y（t）），敵艦位于Q（l，vt）.由于魚雷頭始終對準(zhǔn)敵艦，故此時直線PQ就是魚雷的軌跡曲線弧OP在點P處的切線，即有y=，即=（1—x）y+y······（1）又根據(jù)題意，弧OP的長度的2倍，即dx=2······（2）由（1）（2）消去t，整理得模型：（1—x）y=······（3）六、模型求解解法一：（解析法）特點：有較強(qiáng)的技巧性通過降階法將二階微分方程（3）化為一階微分方程：設(shè)y=p，即dy=p，則y=dp（1—x）dp=由分離變量法有：dp=對等式左邊積分有：dp=ln（p+）+C對等式右邊積分有：=—+Cln（p+）=—+C由題意知y（0）=0，即p（0）=0C=0ln（p+）=—y+=（1—x）通過移項，對兩邊同時平方，可消去得：y==兩邊同時積分得：y=（1—x）—（1—x）+C再由y（0）=0解得魚雷的運動軌跡方程為：y=（1—x）—（1—x）+當(dāng)x=1時，y=，即當(dāng)敵艦航行到點（1，）處時被魚雷擊中，被擊中時間為：t==，若=1，則在t時被擊中解法二：（建立參數(shù)方程求數(shù)值解）特點：求解簡單方便設(shè)t時刻敵艦的坐標(biāo)為（X（t），Y（t）），魚雷的坐標(biāo)為（x（t），y（t））.設(shè)魚雷的速度恒為ω，則（）+（）=ω······（1）由于魚雷頭始終對準(zhǔn)敵艦，故魚雷的速度平行于敵艦與魚雷頭位置的差向量即：（）=λ（），λ>0······（2）消去λ得······（3）因敵艦以速度v沿直線x=1運動，由題意v=0.45，則ω=0.9，X=1，Y=t，因此魚雷的運動軌跡的參數(shù)方程為：用MATLAB實現(xiàn)運算得結(jié)果：y=000.00010.00000.00010.00000.00020.00000.00020.00000.00050.00000.00070.00000.00100.00000.00120.00000.00250.00000.00370.00000.00500.00000.00620.00000.01250.00000.01880.00010.02510.00020.03130.00020.06270.00100.09410.00230.12550.00410.15680.00650.22620.01390.29530.02440.36370.03820.43140.05550.51530.08290.59680.11690.67470.15820.74810.20710.79210.24230.83310.28080.87070.32270.90450.36770.92410.39810.94170.42970.95700.46240.97010.49610.97740.51740.98360.53900.98860.56090.99280.58290.99570.59780.99780.61260.99800.62710.99850.64160.99930.64650.99980.65120.99950.65580.99950.66041.00400.66481.00620.6675由以上的數(shù)據(jù)可知：當(dāng)x=1時，yt=y/0.45=0.667/0.45圖形如下：七、結(jié)果分析、檢驗：此模型將魚雷擊艦問題轉(zhuǎn)化為二圍平面問題，將問題簡單化，利用幾何關(guān)系---曲線的相切性質(zhì)列出常微分方程，利用微分與積分學(xué)求解方程，簡化了計算過程，快速方便的得出所需結(jié)果。模型中也存在不足之處，為了方便模型的建立與問題的求解，我們將魚雷的質(zhì)量忽略不記，故在實際問題中發(fā)射魚雷需水平面以一定的夾角發(fā)射。在魚雷及敵艦運動過程中我們也忽略了水流和氣流向?qū)λ鼈冞\動軌跡的影響。因此所得數(shù)據(jù)與實際數(shù)據(jù)存在一定的誤差。運用MATLAB數(shù)學(xué)軟件避免了復(fù)雜的求解過程，通過程序直接得到了我們所需要的結(jié)果，并且結(jié)果與用解析法求出來的結(jié)果十分接近，再次驗證了結(jié)果的正確性。此模型實用性較強(qiáng)，還可用于其他的追擊問題，如導(dǎo)彈追蹤問題等等。應(yīng)用拓展：我們還可對該問題推廣：如果發(fā)射魚雷時，敵艇立即由儀器覺察，假定敵艇為一高速快艇，他即刻以V的速度與魚雷方向成一夾角逃逸，假設(shè)魚雷的速度為V，問導(dǎo)彈何時何地?fù)糁袛惩?？仿真模型的建立：在這個問題上，就是一步步地模擬導(dǎo)彈追蹤敵艦的實際過程。將該過程分為無數(shù)個細(xì)小的時間微元。將曲線轉(zhuǎn)換為無數(shù)細(xì)小的折線段，再利用我們熟悉的一階數(shù)學(xué)表達(dá)式將魚雷與敵艦的坐標(biāo)表示出來。再利用迭代的思想，依次列出下一時間單元的表達(dá)式，從而推導(dǎo)出其迭代方程，再利用C語言程序求解。設(shè)：逃逸方向與導(dǎo)彈速度方向夾角為θ，如下圖建立坐標(biāo)軸：模型的建立：考慮到艦艇的體積，我們認(rèn)為當(dāng)導(dǎo)彈坐標(biāo)點與艦艇坐標(biāo)點距離小于一米時擊中艦艇。設(shè)：艦艇沿北偏東度航行，魚雷經(jīng)過t時間擊中敵艦，此時不妨設(shè)魚雷坐標(biāo)為（xd，yd），敵艦坐標(biāo)（xj，yj）則由題意：當(dāng)1時，則算擊中。t=時，xd=v，yd=0xj=vcos+1，yj=vsint=2時，=arctanxd=xd+vcos，yd=yd+vsinxj=xj+vcos（+），yj=yj+vsin（+）············由此可得迭代公式xd=v，yd=0xj=vcos+1，yj=vsina=arctanxd=xd+vcosyd=yd+vsinxj=xj+vcos（+）yj=yj+vsin（+）C語言源程序：#include<stdio.h>

#include<math.h>

voidmain()

{

doublex,xd1,xd2,yd1,yd2,v1=0.45,v2=0.9,t,a,q;

intcount=1;

printf("請輸入t和q的值,中間用逗號隔開：");

scanf("%f","%f",&t,&q);

xd1=v1*t;

yd1=0;

xd2=v2*t*cos(q)+1;

yd2=v2*t*sin(q);

x=(yd1-yd2)/(xd1-xd2);

a=atan(x);

do

{

xd1=xd1+v1*t*cos(a);

yd1=yd1+v1*t*sin(a);

xd2=xd2+v2*t*cos(a+q);

yd2=yd2+v2*t*sin(a+q);

x=(yd1-yd2)/(xd1-xd2);

a=atan(x);

count++;

}

while(sqrt((xd1*xd1-xd2*xd2)+(yd1*yd1-yd2*yd2))<=1);

printf("%f",t*count);

}

八、模型的評價與改進(jìn)：該模型真實地模擬了題目的意思，并且通過多種不同的方法求解出非常相近的結(jié)果，同時驗證出了模型的正確性。但是模型本身是基于假設(shè)魚雷的質(zhì)量忽略不記，在魚雷及敵艦運動過程中忽略水流和氣流向?qū)λ鼈冞\動軌跡的影響。所以與實際是有一定的差距的。同時模型中引用了MATLAB軟件來進(jìn)行運算，大大簡化了運算過程，但是其本身是存在誤差的。改進(jìn)：現(xiàn)在我們對其做一些改進(jìn)，在附錄的源程序中，按二分法逐步修改tf，即分別取tf=1，0．5,0．25,…,直到tf=0．21時，得圖2數(shù)據(jù)見附錄中：由圖與數(shù)據(jù)可知改進(jìn)后的模型收斂速度更快。結(jié)論：魚雷需要1.41min才會擊中敵艦九、參考文獻(xiàn)《高等數(shù)學(xué)》第六版上冊，同濟(jì)大學(xué)數(shù)學(xué)系編，高等教育出版社；《MATLAB7.0實用教程》，張圣勤編著，機(jī)械工業(yè)出版社。十、附錄參數(shù)方程解數(shù)值源程序：functiondy=eq2(t,y)dy=zeros(2,1);dy(1)=2*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2);dy(2)=2*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2);t0=0,tf=2;[t,y]=ode45('eq2',[t0tf],[00]);X=1;Y=0:0.01:2;plot(X,Y,'-')holdonplot(y(:,1),y(:,2),'*')改進(jìn)模型的數(shù)據(jù)：000.00010.00000.00010.00000.00020.00000.00020.00000.00050.00000.00070.00000.00100.00000.00120.00000.00250.00000.00370.00000.00500.00000.00620.00000.01250.00000.01880.00010.02510.00020.03130.00020.06270.00100.09410.00230.12550.00410.15680.00650.19160.00980.22640.01390.26100.01880.29560.02440.33000.03090.36420