数值计算实习

1、 数值计算 课程实习报告学 校西南林业大学实习时间2012-12-22专业信息与计算科学班级信息与计算科学2010级分组成员学 号实习题目名称足球比赛中的吊门初射角问题实习条件硬件环境Windows XP软件环境MATLAB 7.1参考文献【1】赵静，但琦，数学建模与数学实验（第2版）。高等教育出版社，2008.1【2】冉启康 张振宇 张立柱，常用数学软件教程。人民邮电出版社，2008.10【3】张德丰，数值分析与应用。国防工业出版社，2007.1【4】郑汉鼎，刁在筠，数学规划M。山东，山东教育出版社，1997.12【5】马正飞，数学计算方法与软件的工程应用。化学工业出版社，2002.12【6

2、】黄有谦，计算方法。北京，高等教育出版社，1981实习内容前言、问题描述、具体理论知识点、具体实例、程序清单、程序实现的核心代码、小组成员分工合作清单。实习过程1、 数值结果2、 出现的问题，解决方法及体会评价1. 写实验报告，占 20 2. 按照教学计划的实验：现场编程序，演示计算结果，占 60 3. 将计算结果用数学软件作图和分析，占 10 4. 与教学内容相关、自由选题或参与教师科研的选题，占 10 2012年 12 月 日足球比赛中的吊门初射角问题 实习内容以及过程前言：吊射是球员在距离球门很远的一种射门方式，一般是看到对方门将站位靠前或者注意力有所分散时做吊射选择，由于足球在空中运行

3、时间较长，可以给守门员反应的时间较长，需要精确的脚法才能成功，因此在职业比赛中出现吊射的情况很罕见，但吊射往往能起到出其不意，打破僵局的作用。吊射的运动轨迹与弹道轨迹相似，可能的影响因素有球与球门的距离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、几组数据分别给出能吊门成功的相应初射角范围，画出最小及最大角度的运动轨迹。要求精度在小数点后至少第3位。(所给数据见计算结果)。· a=6m，b=1m；· a=10m，b=3m；· a=20m，b=5m。 问题分析:1.先考虑最简单情形，即不考虑空气阻力等，此时，球的运动轨迹是抛物线，如果守门员不动，总有合适的角度使吊门成功。2.这不是求一个角度值，而是求一个范围！通常的思路是把问题整理成两个方程求根问题：一个方程是求吊门成功的最小角度，一个方程是求吊门成功的最大角度。3.有可能落地弹入球门，要考虑反弹入门的情况。直观分析:1.最简单情形，抛射体的运动轨迹为抛物线

5、方程如下2.借助于使用方便的数学软件，可直观地看到各种初射角对应的抛射体运动的轨迹图形。 最简情形程序：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; %取整数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)程序抛射体轨迹函数（算法）：fun

7、ction y,t=paosheti1(x,alpha,v,g)% 函数返回抛射体轨迹及第一第一次落地飞行时间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运行结果:1.对于第一组数据，吊门成功的最小角度1.5369（为弧度，下同），对应的时间大约在4.9179秒，最大角度1.53791，对应的时间是5.0689秒；· a=6m，b=1m；最小角度：最大角度2.对于第二组数据，吊门成功的最小角度1.51441，对应的时间大约在4.1403秒，最大角度1.5159，对应的时间大约在4.2526秒；· a=10m，b=3m；最小角度：最大角度：3.对于第三组数据，吊门成功的最小角度1.457，对应的时间大约在4.4033秒，最大角度1.4603，对应的时间大约是4.5343秒· a=20m，b

9、=5m。 最小角度： 最大角度：初步结果分析以及问题的再分析1.遵循由简单到一般的建模原则，先考虑简单情形，得到初步结果，以此为基础，发现问题、分析问题，找到求解思路，并逐步将问题一般化，甚至可以发现逐步一般化的顺序（还是由简到繁，先将哪些方面使之更一般化，如空气阻力、守门员移动）；2.结果有一定的合理性；3.从近似计算角度分析，在允许的精度范围内，如上的“作图观察调整”不失为一种求解方法；注：相比球与守门员及与球门的距离，注意守门员移动的时间，显然守门员有足够的时间移动，因此吊门是不会成功的！ 原因在于将问题假设得过于简单化、理想化了！修改假设应是下一步必须考虑的。小结：最简情形之假设:1.

10、不考虑空气阻力；2.不考虑守门员在球运行过程中的移动；3.球落地是完全弹性的，只考虑仅有一次触地反弹形成的吊门情况；4.只考虑越过守门员头顶的吊门，即出球点与守门员连成一线延伸到球门这样一个直线方向，不考虑从守门员侧面吊门的情况；5.将球看作是数学上的一个点；6.不考虑球的旋转，实际比赛时，旋转是很重要的！7.球的质量为一个单位。有空气阻力的情形（一）仅x方向考虑空气阻力Ø 假设只考虑x方向受空气阻力的影响；Ø 假设空气阻力与速度成正比，比例系数为k=0.4。 Ø 此时，x(t)满足如下的微分方程初值问题问题的解问题的解飞行轨迹：飞行到守门员位置的时间tH：飞行到

11、球门线位置的时间TH：空气阻力的情形（一）（为了考虑时间不足，现在只针对第三组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 空气阻力的情形（一）的结果及分析改进1前面结果有问题，反弹后的角度不应该是alpha了，应该以落地时的情况计算出新反射角。2.修改抛射体函数：将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空气阻力的情形（一）的进一步结果及分析改进Ø 针对第三组数据，计算的最小角度为1.2684，守门员移动时间为

19、2.7836秒，最大角度是1.2807，时间是3.0032秒；最小角度： 最大角度：Ø 结果仍有问题：反弹前后的两波高度一样；Ø 解决的办法是再考虑y方向也有空气阻力。有空气阻力的情形（二）x、y方向均考虑空气阻力Ø 假设x，y两个方向均受空气阻力的影响 ；Ø 假设空气阻力与速度成正比，比例系数为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; %应为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 %私有函数 （为求解t0定义的两个私有函数）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;空气阻力情形（二）的结果分析及改进1.针对第三组数据，计算的最小角度为1.2394，守门员移动时间为2.3843秒，最大角度是1.2482，时间是2.4914秒最小角度：最大角度：2.有必要考虑守门员可以移动的情形。