应用MATLAB建模实例-吊门过程仿真

应用MATLAB建模实例

——数学也是一门技术

数学是一门技术一个例子——足球比赛中的吊门问题谈谈数学建模竞赛及培训数学是一门技术技术的定义《辞海》：泛指根据生产实践经验和自然科学原理而发展成的各种工艺操作方法与技能；除操作技能外，广义的还包括相应的生产工具和其他物质设备，以及生产的工艺过程或作业程序、方法。《科学学辞典》和《科技辞典》：是为社会生产和人类物质文化生活需要服务的，供人类利用和改造自然的物质手段、智能手段和信息手段的总和。数学及其应用的特征是一种智能形态的技术在数学软件的平台上，又表现为一般的物化形态数学建模技术指数学及其应用于解决实际问题的整个过程多样性、合理性、具体问题具体分析艺术性足球比赛中的吊门问题

考虑如下的因素：球与球门的距离为a，守门员与球门的距离为b，球门高h，守门员最大摸高H，球出脚的初速度为v，与水平方向的夹角为alpha（称为初射角）．给定，h=2.44m，H=3.20m，v=30m/s，重力加速度g=10m/s2，针对下列几组数据分别给出能吊门成功的相应初射角范围，要求精度在小数点后第3位。a=6m，b=1m；a=10m，b=3m；a=20m，b=5m。问题分析先考虑最简单情形，即不考虑空气阻力等…，此时，球的运动轨迹是抛物线，如果守门员不动，总有合适的角度使吊门成功。这不是求一个角度值，而是求一个范围！通常的思路是把问题整理成两个方程求根问题：一个方程是求吊门成功的最小角度，一个方程是求吊门成功的最大角度。有可能落地弹入球门，要考虑反弹入门的情况。直观分析最简单情形，抛射体的运动轨迹为抛物线方程如下借助于使用方便的数学软件，可直观地看到各种初射角对应的抛射体运动的轨迹图形。最简情形——程序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之抛射体轨迹函数function[y,t]=paosheti1(x,alpha,v,g)%函数返回抛射体轨迹及第一第一次落地飞行时间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初步结果对于第一组数据，吊门成功的最小角度1.53697（为弧度，下同），对应的时间大约在4.9281秒，最大角度1.53787，对应的时间是5.0627秒；对于第二组数据，吊门成功的最小角度1.51437，对应的时间大约在4.1374秒，最大角度1.51587，对应的时间大约在4.2503秒；对于第三组数据，吊门成功的最小角度1.45718，对应的时间大约在4.4103秒，最大角度1.46022，对应的时间大约是4.531秒。初步结果分析——问题的再分析遵循由简单到一般的建模原则，先考虑简单情形，得到初步结果，以此为基础，发现问题、分析问题，找到求解思路，并逐步将问题一般化，甚至可以发现逐步一般化的顺序（还是由简到繁，先将哪些方面使之更一般化，如空气阻力、守门员移动）；结果有一定的合理性；从近似计算角度分析，在允许的精度范围内，如上的“作图——观察——调整”不失为一种求解方法；相比球与守门员及与球门的距离，注意守门员移动的时间，显然守门员有足够的时间移动，因此吊门是不会成功的！原因在于将问题假设得过于简单化、理想化了！修改假设应是下一步必须考虑的。小结：最简情形之假设不考虑空气阻力；不考虑守门员在球运行过程中的移动；球落地是完全弹性的，只考虑仅有一次触地反弹形成的吊门情况；只考虑越过守门员头顶的吊门，即出球点与守门员连成一线延伸到球门这样一个直线方向，不考虑从守门员侧面吊门的情况；将球看作是数学上的一个点；不考虑球的旋转，实际比赛时，旋转是很重要的！球的质量为一个单位。有空气阻力的情形之一

——仅x方向考虑空气阻力假设只考虑x方向受空气阻力的影响；假设空气阻力与速度成正比，比例系数为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了，应该以落地时的情况计算出新反射角。修改抛射体函数：将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空气阻力的情形之一的结果及分析改进（二）针对第三组数据，计算的最小角度为1.268，守门员移动时间为2.7771秒，最大角度是1.27，时间是2.8101秒；结果仍有问题：反弹前后的两波高度一样；解决的办法是再考虑y方向也有空气阻力。有空气阻力的情形之二

——x、y方向均考虑空气阻力假设x，y两个方向均受空气阻力的影响；假设空气阻力与速度成正比，比例系数为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%私有函数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定义的两个私有函数空气阻力情形之二的结果分析及改进针对第三组数据，计算的最小角度为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;%绘制球移动和守门员移动的动态影像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.