數(shù)值計算實習(xí)

1、 數(shù)值計算 課程實習(xí)報告學(xué) 校西南林業(yè)大學(xué)實習(xí)時間2012-12-22專業(yè)信息與計算科學(xué)班級信息與計算科學(xué)2010級分組成員學(xué) 號實習(xí)題目名稱足球比賽中的吊門初射角問題實習(xí)條件硬件環(huán)境Windows XP軟件環(huán)境MATLAB 7.1參考文獻【1】趙靜，但琦，數(shù)學(xué)建模與數(shù)學(xué)實驗（第2版）。高等教育出版社，2008.1【2】冉啟康 張振宇 張立柱，常用數(shù)學(xué)軟件教程。人民郵電出版社，2008.10【3】張德豐，數(shù)值分析與應(yīng)用。國防工業(yè)出版社，2007.1【4】鄭漢鼎，刁在筠，數(shù)學(xué)規(guī)劃M。山東，山東教育出版社，1997.12【5】馬正飛，數(shù)學(xué)計算方法與軟件的工程應(yīng)用。化學(xué)工業(yè)出版社，2002.12【6

2、】黃有謙，計算方法。北京，高等教育出版社，1981實習(xí)內(nèi)容前言、問題描述、具體理論知識點、具體實例、程序清單、程序?qū)崿F(xiàn)的核心代碼、小組成員分工合作清單。實習(xí)過程1、 數(shù)值結(jié)果2、 出現(xiàn)的問題，解決方法及體會評價1. 寫實驗報告，占 20 2. 按照教學(xué)計劃的實驗：現(xiàn)場編程序，演示計算結(jié)果，占 60 3. 將計算結(jié)果用數(shù)學(xué)軟件作圖和分析，占 10 4. 與教學(xué)內(nèi)容相關(guān)、自由選題或參與教師科研的選題，占 10 2012年 12 月 日足球比賽中的吊門初射角問題 實習(xí)內(nèi)容以及過程前言：吊射是球員在距離球門很遠的一種射門方式，一般是看到對方門將站位靠前或者注意力有所分散時做吊射選擇，由于足球在空中運行

3、時間較長，可以給守門員反應(yīng)的時間較長，需要精確的腳法才能成功，因此在職業(yè)比賽中出現(xiàn)吊射的情況很罕見，但吊射往往能起到出其不意，打破僵局的作用。吊射的運動軌跡與彈道軌跡相似，可能的影響因素有球與球門的距離a，守門員與球門的距離b，球門高h，守門員最大接球高度H，球在空中飛行時間t，球出腳的初速度v，與水平方向的夾角alpha（初射角）。問題的提出與分析考慮如下的因素：球與球門的距離為a，守門員與球門的距離為b，球門高h，守門員最大摸高H，球出腳的初速度為v，與水平方向的夾角為alpha（稱為初射角）給定，h=2.44m，H=3.20m，v=30m/s，重力加速度g=10m/s2提出問題：針對下列

4、幾組數(shù)據(jù)分別給出能吊門成功的相應(yīng)初射角范圍，畫出最小及最大角度的運動軌跡。要求精度在小數(shù)點后至少第3位。(所給數(shù)據(jù)見計算結(jié)果)。· a=6m，b=1m；· a=10m，b=3m；· a=20m，b=5m。 問題分析:1.先考慮最簡單情形，即不考慮空氣阻力等，此時，球的運動軌跡是拋物線，如果守門員不動，總有合適的角度使吊門成功。2.這不是求一個角度值，而是求一個范圍！通常的思路是把問題整理成兩個方程求根問題：一個方程是求吊門成功的最小角度，一個方程是求吊門成功的最大角度。3.有可能落地彈入球門，要考慮反彈入門的情況。直觀分析:1.最簡單情形，拋射體的運動軌跡為拋物線

5、方程如下2.借助于使用方便的數(shù)學(xué)軟件，可直觀地看到各種初射角對應(yīng)的拋射體運動的軌跡圖形。 最簡情形程序：v=30;g=10; h=2.44;H=3.2;a=6;b=1; l=a-b;L=a*1.01; %1.01表示進入門里x=0:0.01:L; %對x采樣i1=floor(a/0.01)+1; i2=floor(a-b)/0.01)+1; %取整數(shù)alpha=1.5368:0.00001:1.538; n=length(alpha);%弧度for i=1:n; y,tfinal=paosheti1(x,alpha(i),v,g); tH=l/(v*cos(alpha(i); %從射門到球到守

6、門員位置的時間 if y(i1)>=h | y(i2)<=H; success(i)=0; else success(i)=1; endplot(l,H,'r+',a,h,'r+'), hold on, plot(x,y), grid, hold off title('足球比賽中的吊門 ','初射角=',num2str(alpha(i),6) ,. ' 守門員的移動時間=',num2str(tH),pauseendfigure(2)plot(alpha,success)程序拋射體軌跡函數(shù)（算法）：fun

7、ction y,t=paosheti1(x,alpha,v,g)% 函數(shù)返回拋射體軌跡及第一第一次落地飛行時間y=x*tan(alpha)-x.2*g/(2*v2*(cos(alpha)2);% 不考慮是直接進門還是落地后進門的軌跡t=2*v*sin(alpha)/g; % 第一次落地前最大飛行時間xmax=v*cos(alpha)*t; %第一次落地前最大飛行X向距離n=length(x); for i=1:n if y(i)<0 %判斷是落地后進門 xx=x(i)-xmax; y(i)=xx*tan(alpha)-xx.2*g/(2*v2*(cos(alpha)2); %計算落地后進

8、門的后半段軌跡 endend運行結(jié)果:1.對于第一組數(shù)據(jù)，吊門成功的最小角度1.5369（為弧度，下同），對應(yīng)的時間大約在4.9179秒，最大角度1.53791，對應(yīng)的時間是5.0689秒；· a=6m，b=1m；最小角度：最大角度2.對于第二組數(shù)據(jù)，吊門成功的最小角度1.51441，對應(yīng)的時間大約在4.1403秒，最大角度1.5159，對應(yīng)的時間大約在4.2526秒；· a=10m，b=3m；最小角度：最大角度：3.對于第三組數(shù)據(jù)，吊門成功的最小角度1.457，對應(yīng)的時間大約在4.4033秒，最大角度1.4603，對應(yīng)的時間大約是4.5343秒· a=20m，b

9、=5m。 最小角度： 最大角度：初步結(jié)果分析以及問題的再分析1.遵循由簡單到一般的建模原則，先考慮簡單情形，得到初步結(jié)果，以此為基礎(chǔ)，發(fā)現(xiàn)問題、分析問題，找到求解思路，并逐步將問題一般化，甚至可以發(fā)現(xiàn)逐步一般化的順序（還是由簡到繁，先將哪些方面使之更一般化，如空氣阻力、守門員移動）；2.結(jié)果有一定的合理性；3.從近似計算角度分析，在允許的精度范圍內(nèi)，如上的“作圖觀察調(diào)整”不失為一種求解方法；注：相比球與守門員及與球門的距離，注意守門員移動的時間，顯然守門員有足夠的時間移動，因此吊門是不會成功的！ 原因在于將問題假設(shè)得過于簡單化、理想化了！修改假設(shè)應(yīng)是下一步必須考慮的。小結(jié)：最簡情形之假設(shè):1.

10、不考慮空氣阻力；2.不考慮守門員在球運行過程中的移動；3.球落地是完全彈性的，只考慮僅有一次觸地反彈形成的吊門情況；4.只考慮越過守門員頭頂?shù)牡蹰T，即出球點與守門員連成一線延伸到球門這樣一個直線方向，不考慮從守門員側(cè)面吊門的情況；5.將球看作是數(shù)學(xué)上的一個點；6.不考慮球的旋轉(zhuǎn)，實際比賽時，旋轉(zhuǎn)是很重要的！7.球的質(zhì)量為一個單位。有空氣阻力的情形（一）僅x方向考慮空氣阻力Ø 假設(shè)只考慮x方向受空氣阻力的影響；Ø 假設(shè)空氣阻力與速度成正比，比例系數(shù)為k=0.4。 Ø 此時，x(t)滿足如下的微分方程初值問題問題的解問題的解飛行軌跡：飛行到守門員位置的時間tH：飛行到

11、球門線位置的時間TH：空氣阻力的情形（一）（為了考慮時間不足，現(xiàn)在只針對第三組a=20;b=5時的情況，下同）程序v=30;k=0.4;g=10;h=2.44;H=3.2;a=20;b=5;l=a-b;L=a*1.1;for alpha=1.2:0.001:1.3 %1.5425%pi/2-eps Th=-log(1-a*k/(v*cos(alpha)/k; T=Th*1.2; t=0:0.01:T; x,y=paosheti2(t,alpha,v,k,g); TH=-log(1-l*k/(v*cos(alpha)/k; plot(l,H,'r+',a,h,'r+

12、9;),hold on, plot(x,y),grid, hold off title('足球比賽中的吊門 ','初射角=',num2str(alpha,6) ,. ' 守門員的移動時間=',num2str(TH),pauseend 空氣阻力的情形（一）算法2：function x,y=paosheti2(t,alpha,v,k,g)x=v*cos(alpha)*(1-exp(-k*t)/k;y=v*sin(alpha)*t-g*t.2/2;n=length(t);t0=2*v*sin(alpha)/g;%the time when the ba

13、ll down to groundxt0=v*cos(alpha)*(1-exp(-k*t0)/k;vxt0=v*cos(alpha)*exp(-k*t0);vyt0=v*sin(alpha);vt0=sqrt(vxt02+vyt02);for i=1:n if t(i)>t0 tt=t(i)-t0; x(i)=xt0+vt0*cos(alpha)*(1-exp(-k*tt)/k; y(i)=vt0*sin(alpha)*tt-g*tt2/2; endend function x,y=paosheti2(t,alpha,v,k,g)%返回拋物體軌跡（考慮空氣阻力）x=v*cos(alph

14、a)*(1-exp(-k*t)/k;y=v*sin(alpha)*t-g*t.2/2;n=length(t);t0=2*v*sin(alpha)/g; % 第一次落地前最大飛行時間xt0=v*cos(alpha)*(1-exp(-k*t0)/k; % 第一次落地前X向最大飛行距離vxt0=v*cos(alpha)*exp(-k*t0); % 第一次落地X向飛行速度vyt0=v*sin(alpha); % 第一次落地Y向飛行速度vt0=sqrt(vxt02+vyt02); % 第一次落地飛行速度大小for i=1:n if t(i)>t0 % 判斷是落地后進門 tt=t(i)-t0; %

15、計算落地后進門的后半段軌跡 x(i)=xt0+vt0*cos(alpha)*(1-exp(-k*tt)/k; y(i)=vt0*sin(alpha)*tt-g*tt2/2; endend 空氣阻力的情形（一）的結(jié)果及分析改進1前面結(jié)果有問題，反彈后的角度不應(yīng)該是alpha了，應(yīng)該以落地時的情況計算出新反射角。2.修改拋射體函數(shù)：將paosheti2(t,alpha,v,k,g)，換成paosheti22(t,alpha,v,k,g)。 算法2修改如下得到算法22function x,y=paosheti22(t,alpha,v,k,g)x=v*cos(alpha)*(1-exp(-k*t)/k

16、;y=v*sin(alpha)*t-g*t.2/2;n=length(t);t0=2*v*sin(alpha)/g;%the time when the ball down to the groundxt0=v*cos(alpha)*(1-exp(-k*t0)/k;vxt0=v*cos(alpha)*exp(-k*t0);vyt0=v*sin(alpha);vt0=sqrt(vxt02+vyt02);alpha1=atan(vyt0/vxt0);for i=1:n if t(i)>t0 tt=t(i)-t0; x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt)/

17、k; y(i)=vt0*sin(alpha1)*tt-g*tt2/2; endend function x,y=paosheti22(t,alpha,v,k,g)%返回拋物體軌跡（考慮空氣阻力、反彈角度）x=v*cos(alpha)*(1-exp(-k*t)/k;y=v*sin(alpha)*t-g*t.2/2;n=length(t);t0=2*v*sin(alpha)/g; % 第一次落地前最大飛行時間xt0=v*cos(alpha)*(1-exp(-k*t0)/k; % 第一次落地前X向最大飛行距離vxt0=v*cos(alpha)*exp(-k*t0); % 第一次落地X向飛行速度vyt

18、0=v*sin(alpha); % 第一次落地Y向飛行速度vt0=sqrt(vxt02+vyt02); % 第一次落地飛行速度大小alpha1=atan(vyt0/vxt0); % 計算反彈角度for i=1:n if t(i)>t0 % 判斷是落地后進門 tt=t(i)-t0; % 計算落地后進門的后半段軌跡 x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt)/k; y(i)=vt0*sin(alpha1)*tt-g*tt2/2; endend空氣阻力的情形（一）的進一步結(jié)果及分析改進Ø 針對第三組數(shù)據(jù)，計算的最小角度為1.2684，守門員移動時間為

19、2.7836秒，最大角度是1.2807，時間是3.0032秒；最小角度： 最大角度：Ø 結(jié)果仍有問題：反彈前后的兩波高度一樣；Ø 解決的辦法是再考慮y方向也有空氣阻力。有空氣阻力的情形（二）x、y方向均考慮空氣阻力Ø 假設(shè)x，y兩個方向均受空氣阻力的影響 ；Ø 假設(shè)空氣阻力與速度成正比，比例系數(shù)為k=0.4。 Ø 此時，x(t)仍滿足同上的常微分方程初值問題y(t)滿足如下的常微分方程初值問題 問題的解 有空氣阻力的情形（二）程序：v=30;k=0.4;g=10;h=2.44;H=3.2;a=20;b=5;l=a-b;L=a*1.1;for a

20、lpha=1.2:0.001:1.3 %1.5425%pi/2-eps Th=-log(1-a*k/(v*cos(alpha)/k; T=Th*1.2; t=0:0.01:T; x,y=paosheti3(t,alpha,v,k,g); TH=-log(1-l*k/(v*cos(alpha)/k; plot(l,H,'r+',a,h,'r+'),hold on, plot(x,y),grid, hold off title('足球比賽中的吊門 ','初射角=',num2str(alpha,6) ,. ' 守門員的移動時間=

21、',num2str(TH),pauseend 有空氣阻力的情形（二）算法：function x,y=paosheti3(t,alpha,v,k,g) %返回拋物體軌跡（考慮空氣阻力）x=v*cos(alpha)*(1-exp(-k*t)/k;y=(v*sin(alpha)+g/k)*(1-exp(-k*t)/k-g*t/k;n=length(t); % 第一次落地前最大飛行時間t00=2.; %在2附近尋求零點tt0(1)=t00;tb=1; ii=1;while(abs(tb)>1e-5) %收斂條件 tt0(ii+1)=tt0(ii)-paoshetiy(tt0(ii),al

22、pha,v,k,g)/dpaoshetiy(tt0(ii),alpha,v,k,g); tb=tt0(ii+1)-tt0(ii); ii=ii+1; if(ii>20)error('numb. of iter. is 30 times'); endend t0=tt0(ii); % 收斂到零點， 第一次落地前最大飛行時間y0=(v*sin(alpha)+g/k)*(1-exp(-k*t0)/k-g*t0/k; %應(yīng)為0或接近0xt0=v*cos(alpha)*(1-exp(-k*t0)/k; % 第一次落地前X向最大飛行距離vxt0=v*cos(alpha)*exp(-k

23、*t0); % 第一次落地X向飛行速度vyt0=(v*sin(alpha)+g/k)*exp(-k*t0)-g/k; % 第一次落地Y向飛行速度vt0=sqrt(vxt02+vyt02); % 第一次落地飛行速度大小alpha1=atan(abs(vyt0/vxt0); % 計算反彈角度for i=1:n if t(i)>t0 % 判斷是落地后進門 tt=t(i)-t0; % 計算落地后進門的后半段軌跡 x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt)/k; y(i)=(vt0*sin(alpha1)+g/k)*(1-exp(-k*tt)/k-g*tt/k; endend %私有函數(shù) （為求解t0定義的兩個私有函數(shù)）function y=paoshetiy(t,alpha,v,k,g)y=(v*sin(alpha)+g/k)*(1-exp(-k*t)/k-g*t/k;function dy=dpaoshetiy(t,alpha,v,k,g)dy=(v*sin(alpha)+g/k)*exp(-k*t)-g/k;空氣阻力情形（二）的結(jié)果分析及改進1.針對第三組數(shù)據(jù)，計算的最小角度為1.2394，守門員移動時間為2.3843秒，最大角度是1.2482，時間是2.4914秒最小角度：最大角度：2.有必要考慮守門員可以移動的情形。