應用MATLAB建模實例-吊門過程仿真

應用MATLAB建模實例

——數(shù)學也是一門技術

數(shù)學是一門技術一個例子——足球比賽中的吊門問題談談數(shù)學建模競賽及培訓數(shù)學是一門技術技術的定義《辭?！罚悍褐父鶕?jù)生產實踐經驗和自然科學原理而發(fā)展成的各種工藝操作方法與技能；除操作技能外，廣義的還包括相應的生產工具和其他物質設備，以及生產的工藝過程或作業(yè)程序、方法。《科學學辭典》和《科技辭典》：是為社會生產和人類物質文化生活需要服務的，供人類利用和改造自然的物質手段、智能手段和信息手段的總和。數(shù)學及其應用的特征是一種智能形態(tài)的技術在數(shù)學軟件的平臺上，又表現(xiàn)為一般的物化形態(tài)數(shù)學建模技術指數(shù)學及其應用于解決實際問題的整個過程多樣性、合理性、具體問題具體分析藝術性足球比賽中的吊門問題

考慮如下的因素：球與球門的距離為a，守門員與球門的距離為b，球門高h，守門員最大摸高H，球出腳的初速度為v，與水平方向的夾角為alpha（稱為初射角）．給定，h=2.44m，H=3.20m，v=30m/s，重力加速度g=10m/s2，針對下列幾組數(shù)據(jù)分別給出能吊門成功的相應初射角范圍，要求精度在小數(shù)點后第3位。a=6m，b=1m；a=10m，b=3m；a=20m，b=5m。問題分析先考慮最簡單情形，即不考慮空氣阻力等…，此時，球的運動軌跡是拋物線，如果守門員不動，總有合適的角度使吊門成功。這不是求一個角度值，而是求一個范圍！通常的思路是把問題整理成兩個方程求根問題：一個方程是求吊門成功的最小角度，一個方程是求吊門成功的最大角度。有可能落地彈入球門，要考慮反彈入門的情況。直觀分析最簡單情形，拋射體的運動軌跡為拋物線方程如下借助于使用方便的數(shù)學軟件，可直觀地看到各種初射角對應的拋射體運動的軌跡圖形。最簡情形——程序1-1v=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;alpha=1.5368:0.00001:1.538;n=length(alpha);fori=1:n; [y,tfinal]=paosheti1(x,alpha(i),v,g); tH=l/(v*cos(alpha(i)));%從射門到球到守門員位置的時間 ify(i1)>=h|y(i2)<=H; success(i)=0; else success(i)=1; end plot(l,H,'r+',a,h,'r+'),holdon, plot(x,y),grid,holdoff title(['足球比賽中的吊門','初射角=',num2str(alpha(i),6),...'守門員的移動時間=',num2str(tH)]),pauseendfigure(2)plot(alpha,success)程序1-1之拋射體軌跡函數(shù)function[y,t]=paosheti1(x,alpha,v,g)%函數(shù)返回拋射體軌跡及第一第一次落地飛行時間y=x*tan(alpha)-x.^2*g/(2*v^2*(cos(alpha))^2);%不考慮是直接進門還是落地后進門的軌跡t=2*v*sin(alpha)/g;%第一次落地前最大飛行時間xmax=v*cos(alpha)*t;%第一次落地前最大飛行X向距離n=length(x);fori=1:nify(i)<0%判斷是落地后進門xx=x(i)-xmax;y(i)=xx*tan(alpha)-xx.^2*g/(2*v^2*(cos(alpha))^2);

%計算落地后進門的后半段軌跡endend初步結果對于第一組數(shù)據(jù)，吊門成功的最小角度1.53697（為弧度，下同），對應的時間大約在4.9281秒，最大角度1.53787，對應的時間是5.0627秒；對于第二組數(shù)據(jù)，吊門成功的最小角度1.51437，對應的時間大約在4.1374秒，最大角度1.51587，對應的時間大約在4.2503秒；對于第三組數(shù)據(jù)，吊門成功的最小角度1.45718，對應的時間大約在4.4103秒，最大角度1.46022，對應的時間大約是4.531秒。初步結果分析——問題的再分析遵循由簡單到一般的建模原則，先考慮簡單情形，得到初步結果，以此為基礎，發(fā)現(xiàn)問題、分析問題，找到求解思路，并逐步將問題一般化，甚至可以發(fā)現(xiàn)逐步一般化的順序（還是由簡到繁，先將哪些方面使之更一般化，如空氣阻力、守門員移動）；結果有一定的合理性；從近似計算角度分析，在允許的精度范圍內，如上的“作圖——觀察——調整”不失為一種求解方法；相比球與守門員及與球門的距離，注意守門員移動的時間，顯然守門員有足夠的時間移動，因此吊門是不會成功的！原因在于將問題假設得過于簡單化、理想化了！修改假設應是下一步必須考慮的。小結：最簡情形之假設不考慮空氣阻力；不考慮守門員在球運行過程中的移動；球落地是完全彈性的，只考慮僅有一次觸地反彈形成的吊門情況；只考慮越過守門員頭頂?shù)牡蹰T，即出球點與守門員連成一線延伸到球門這樣一個直線方向，不考慮從守門員側面吊門的情況；將球看作是數(shù)學上的一個點；不考慮球的旋轉，實際比賽時，旋轉是很重要的！球的質量為一個單位。有空氣阻力的情形之一

——僅x方向考慮空氣阻力假設只考慮x方向受空氣阻力的影響；假設空氣阻力與速度成正比，比例系數(shù)為k=0.4。此時，x(t)滿足如下的微分方程初值問題問題的解飛行軌跡飛行到守門員位置的時間tH：問題的解飛行軌跡飛行到球門線位置的時間TH：空氣阻力的情形之一——程序2-1-1clc,clearallv=30;g=10;h=2.44;H=3.2;a=20;b=5;k=0.4;l=a-b;L=a*1.01;%1.01表示進入門里x=0:0.01:L;%對x采樣alpha=1.2:0.001:1.3;n=length(alpha);figure(1)fori=1:n; Th=-log(1-a*k/(v*cos(alpha(i))))/k;%計算球從射門到進球門的時間范圍T=Th*1.02;%球過球門線t=0:0.01:T;%對時間采樣[x,y]=paosheti2(t,alpha(i),v,k,g);%計算求得飛行軌跡TH=-log(1-l*k/(v*cos(alpha(i))))/k;%計算從射門到球員位置的時間i1=floor(Th/0.01)+1;%計算球過球門線時刻的采樣點序號i2=floor(TH/0.01)+1;%計算球過守門員時刻的采樣點序號 ify(i1)>=h|y(i2)<=H;%判定是否射門成功 success(i)=0; else success(i)=1; endplot(l,H,'r+',a,h,'r+'),holdon,plot(x,y),grid,holdofftitle(['足球比賽中的吊門','初射角=',num2str(alpha(i),6),...'守門員的移動時間=',num2str(TH),'成功=',num2str(success(i))]),pauseendfigure(2)plot(alpha,success);空氣阻力的情形之一——程序2-1-2function[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;%第一次落地前最大飛行時間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(vxt0^2+vyt0^2);%第一次落地飛行速度大小fori=1:nift(i)>t0%判斷是落地后進門tt=t(i)-t0;%計算落地后進門的后半段軌跡x(i)=xt0+vt0*cos(alpha)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha)*tt-g*tt^2/2;endend空氣阻力的情形之一的結果及分析改進（一）前面結果有問題，反彈后的角度不應該是alpha了，應該以落地時的情況計算出新反射角。修改拋射體函數(shù)：將paosheti2(t,alpha,v,k,g)，換成paosheti22(t,alpha,v,k,g)。程序2-1-2修改如下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向飛行速度vyt0=v*sin(alpha);%第一次落地Y向飛行速度vt0=sqrt(vxt0^2+vyt0^2);%第一次落地飛行速度大小alpha1=atan(vyt0/vxt0);

%計算反彈角度fori=1:nift(i)>t0%判斷是落地后進門tt=t(i)-t0;%計算落地后進門的后半段軌跡x(i)=xt0+vt0*cos(alpha1)*(1-exp(-k*tt))/k;y(i)=vt0*sin(alpha1)*tt-g*tt^2/2;endend空氣阻力的情形之一的結果及分析改進（二）針對第三組數(shù)據(jù)，計算的最小角度為1.268，守門員移動時間為2.7771秒，最大角度是1.27，時間是2.8101秒；結果仍有問題：反彈前后的兩波高度一樣；解決的辦法是再考慮y方向也有空氣阻力。有空氣阻力的情形之二

——x、y方向均考慮空氣阻力假設x，y兩個方向均受空氣阻力的影響；假設空氣阻力與速度成正比，比例系數(shù)為k=0.4。此時，x(t)仍滿足同上的常微分方程初值問題y(t)滿足如下的常微分方程初值問題問題的解有空氣阻力的情形之二——程序2-2-1v=30;k=0.4;g=10;h=2.44;H=3.2;a=20;b=5;l=a-b;L=a*1.1;foralpha=1.2:0.001:1.3%1.5425%pi/2-epsTh=-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+'),holdon,plot(x,y),grid,holdofftitle(['足球比賽中的吊門','初射角=',num2str(alpha,6),...'守門員的移動時間=',num2str(TH)]),pauseend有空氣阻力的情形之二——程序2-2-2function[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),alpha,v,k,g)/dpaoshetiy(tt0(ii),alpha,v,k,g);tb=tt0(ii+1)-tt0(ii);ii=ii+1;if(ii>20)error('numb.ofiter.is30times');endendt0=tt0(ii);%收斂到零點，第一次落地前最大飛行時間y0=(v*sin(alpha)+g/k)*(1-exp(-k*t0))/k-g*t0/k;%應為0或接近0xt0=v*cos(alpha)*(1-exp(-k*t0))/k;%第一次落地前X向最大飛行距離vxt0=v*cos(alpha)*exp(-k*t0);%第一次落地X向飛行速度vyt0=(v*sin(alpha)+g/k)*exp(-k*t0)-g/k;%第一次落地Y向飛行速度vt0=sqrt(vxt0^2+vyt0^2);%第一次落地飛行速度大小alpha1=atan(abs(vyt0/vxt0));%計算反彈角度fori=1:nift(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ù)functiony=paoshetiy(t,alpha,v,k,g)y=(v*sin(alpha)+g/k)*(1-exp(-k*t))/k-g*t/k;functiondy=dpaoshetiy(t,alpha,v,k,g)dy=(v*sin(alpha)+g/k)*exp(-k*t)-g/k;為求解t0定義的兩個私有函數(shù)空氣阻力情形之二的結果分析及改進針對第三組數(shù)據(jù)，計算的最小角度為1.239，守門員移動時間為2.3797秒，最大角度是1.248，時間是2.4889秒；有必要考慮守門員可以移動的情形。守門員可以移動的情形假設守門員只沿球運行的方向移動；球一出腳守門員即判斷好并移動；守門員的移動速度記為u，其大小是吊門能否成功的另一關鍵因素。問題分析球的運動軌跡與上一種情況完全一致，程序2-2-2不必改。單獨加上守門員的移動顯示，二者疊加，可得守門員可移動情形下吊門成功與否的直觀顯示。守門員可以移動的情形之程序3-1-1clc,clearallv=30;g=10;h=2.44;H=3.2;a=20;k=0.4;b0=5;u0=0.5;l0=a-b0;L=a*1.01;%1.01表示進入門里alpha=1.2:0.01:1.3;n=length(alpha);figure(1)fori=1:n; Th=-log(1-a*k/(v*cos(alpha(i))))/k;%計算球從射門到進球門的時間范圍T=Th*1.02;%球過球門線t=0:0.1:T;%對時間采樣[x,y]=paosheti3(t,alpha(i),v,k,g);%計算求得飛行軌跡nn=length(t);forj=1:nn;%繪制球移動和守門員移動的動態(tài)影像xx=x(1:j);yy=y(1:j);ll0=l0+u0*t(j);ifll0>al(j)=a;elsel(j)=ll0;endTH=-log(1-l(j)*k/(v*cos(alpha(i))))/k;%計算球從射門到球員移動的當前位置的時間plot(xx,yy,'bo',l(j),H,'b+',a,h,'ro'),grid,pauseendfori=1:n;

…….%接上部分，分析繪制球飛行軌跡和守門員運動軌跡i1=floor(Th/0.1)+1;%計算球過球門線時刻的采樣點序號i2=floor(TH/0.