2019年认证杯SPSSPRO杯数学建模A题(第二阶段)好风凭借力,送我上青云全过程文档及程序

2019年认证杯SPSSPRO杯数学建模

纸飞机飞行最佳投掷角度及最远距离研究

A题 好风凭借力,送我上青云

原题再现：

  纸飞机有许多种折法。世界上有若干具有一定影响力的纸飞机比赛，通常的参赛规定是使用一张特定规格的纸，例如 A4 大小的纸张，折成一架纸飞机。大多数比赛都不允许使用剪刀和胶水，有的比赛中可以少量使用胶带以调整纸飞机的重心。折好后在无风的室内环境投掷，并测量其留空时间和飞行距离等成绩。在 2012 年，一种叫做 Suzanne 的折法曾创下飞行距离的吉尼斯世界纪录，飞行了超过 69 米。Suzanne 的折法参见：https://v.qq.com/x/page/o0636km20dz.html
  第二阶段问题： 我们依然假定纸飞机 (以 Suzanne 为例) 的翼面是精确的平面。假设我们可以使用少量胶带粘贴在机身上，用以调整重心的前后位置。请建立合理的数学模型，计算出为了使飞行距离尽可能远，飞机重心的最佳位置以及投掷出手的最佳角度。如果我们希望留空时间尽可能长，重心位置和出手角度又该如何取值呢？

整体求解过程概述(摘要)

  气流对机翼的冲击作用能够平衡机身的重力这是纸飞机的特点，在不考虑机翼弯曲的情况下，气流对机翼的作用不存在压强差，仅剩气流对机翼的冲击。由动量定理可知气流对机翼的冲击作用与机翼相对空气的流速有关系，由于机翼形状特殊，且难以知道气流冲击机翼之后的速度分布，所以很难算出冲击力，不过一般与流速的指数次成正比，文中把这个指数称为气流指数，是气流与机翼复杂综合作用的结果，一般可以从实验得到，本文还得到水平方向的气流指数和竖直方向的气流指数并不相同，可能受扰流及机身机构的影响。
  对问题的分析：要在竖直方向上获得飞行高度争取飞行时间，与水平方向上速度不变是矛盾的，因此必然存在一个合适的仰角使得飞行距离最远，这正是要求解的问题，为了得到飞行距离与仰角的关系，以 Suzanne 纸飞机抛出时在地面上的投影为坐标原点，纸飞机飞行方向为横坐标，建立平面直角坐标系，这也就假设了在平面内考虑该问题。将 Suzanne 纸飞机的运动在第一象限内正交分解，由牛顿运动定律得到两个方向上的运动学方程，由于上升阶段与下降阶段受力不同，因此需要分开讨论该问题，最终得到的运动学方程是一个分段函数，由分段函数的连续性可以得到水平距离与仰角的函数从而解决该问题。
  对数据的处理：计算结果表明这是一个非常复杂的函数，不能用一般的方法得到解析解，因此需要将角度离散得到对应水平距离的一组数值解，最后拟合离散角度与对应的水平距离数值得到函数图像，从图中找到最远距离对应的角度，本文的结果约为 42.30。
  本文的问题：计算中发现竖直下落的速度太大，也就是冲击力对竖直速度的平方成正比的假设是不对的，同理水平速度过小，均表示气流指数需要重新调整，这也是本文的最远距离（57m）小于 Suzanne 实验距离的原因，同时也是后面需要修正的地方，通过时间切片重构图片，可以将视频中飞行轨迹数据提取出来，与计算轨迹进行比对，是修身气流指数的方法。本文没有从机身的结构去分析气流对机体的作用，显然没有直面问题，这是后续要研究的地方。

问题分析：

  气流对机翼的冲击作用能够平衡机身的重力这是纸飞机的特点，在不考虑机翼弯曲的情况下，气流对机翼的作用不存在压强差，仅剩气流对机翼的冲击。由动量定理可知气流对机翼的冲击作用与机翼相对空气的流速有关系，由于机翼形状特殊，且难以知道气流冲击机翼之后的速度分布，所以很难算出冲击力，不过一般与流速的指数次成正比，文中把这个指数称为气流指数，是气流与机翼复杂综合作用的结果，一般可以从实验得到，本文还得到水平方向的气流指数和竖直方向的气流指数并不相同，可能受扰流及机身机构的影响。
  对问题的分析：要在竖直方向上获得飞行高度争取飞行时间，与水平方向上速度不变是矛盾的，因此必然存在一个合适的仰角使得飞行距离最远，这正是要求解的问题，为了得到飞行距离与仰角的关系，以 Suzanne 纸飞机抛出时在地面上的投影为坐标原点，纸飞机飞行方向为横坐标，建立平面直角坐标系，这也就假设了在平面内考虑该问题。将 Suzanne 纸飞机的运动在第一象限内正交分解，由牛顿运动定律得到两个方向上的运动学方程，由于上升阶段与下降阶段受力不同，因此需要分开讨论该问题，最终得到的运动学方程是一个分段函数，由分段函数的连续性可以得到水平距离与仰角的函数从而解决该问题。

模型假设：

  1. 上升过程中飞行姿势（仰角）没有发生改变。
  2. 在整个飞行过程中机翼始终没有弯曲变形，始终在同一数值平面内飞行。
  3. 飞行姿势在最高点处对时间不可导。
  4. 气流对机翼冲击后的速度为零。

论文缩略图：

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档not free)

%质心坐标计算
xa1=(int(0.75*x^2,0,1.7)+int((0.24*x+2.79)*x,1.7,2)-int((-6.35*x+12.8
3)*x,1.7,2)+int((0.24*x+2.79)*x,2,9.8)+int((-3.71*x+41.5)*x,9.8,11.2)
)/(int(0.75*x,0,1.7)+int(0.24*x+2.79,1.7,2)-int(-6.85*x+12.83,1.7,2)+
int(0.24*x+2.79,2,9.8)+int(-3.71*x+41.5,9.8,11.2))
ya1=(int(((x-41.5)/-3.71-x/0.75)*x,0,1.5)+int(((x-41.5)/-3.71-(x-12.8
3)/-5.6)*x,1.5,3.2)+int(((x-41.5)/-3.71-(x-2.79)/0.24)*x,3.2,5.2))/(i
nt(0.75*x,0,1.7)+int(0.24*x+2.79,1.7,2)-int(-6.85*x+12.83,1.7,2)+int(
0.24*x+2.79,2,9.8)+int(-3.71*x+41.5,9.8,11.2))
xa2=(int(0.81*x^2,0,1.7)+int(-2.77*x^2+6.01*x,1.7,2.2))/(int(0.81*x,0
,1.7)+int(-2.77*x+6.01,1.7,2.2))
ya2=(int((x-7.26)/-3.3-(x-1.25)/-0.57,0,1.25)+int((x-7.26)/-3.3-(x-1.
25)/0.24,1.25,1.67))/(int(0.81*x,0,1.7)+int(-2.77*x+6.01,1.7,2.2))
xa3=(int(2.77*x^2,0,0.5)+int(1.57*x^2+0.6*x,0.5,1.4)+int(-0.62*x^2+3.
67*x,1.4,1.5)+int(-1.37*x^2+4.77*x,1.5,3.5))/(int(2.77*x,0,0.5)+int(1
.57*x+0.6,0.5,1.4)+int(-0.62*x+3.67,1.4,1.5)+int(-1.37*x+4.77,1.5,3.5
))
ya3=(int(((x+1.1)/0.75-(x-0.6)/-1.2)*x,0,0.6)+int(((x+1.1)/0.75-(x-0.
6)/1.57)*x,0.6,1.5)+int(((x-3.67)/-0.62-(x-0.6)/1.57)*x,1.5,2.8))/(in
t(2.77*x,0,0.5)+int(1.57*x+0.6,0.5,1.4)+int(-0.62*x+3.67,1.4,1.5)+int
(-1.37*x+4.77,1.5,3.5))
xa4=(int(1.42*x^2,0,1.3)+int(-1.42*x^2+3.7*x,1.3,2.2))/(int(1.42*x,0,
1.3)+int(-1.42*x+3.7,1.3,2.2))
ya4=(int(((x+2.2)/1.69-(x-1.5)/-1.15)*x,0,1.5)+int(((x+2.2)/1.69-(x-1
.5)/0.27)*x,1.5,2.6))/(int(1.42*x,0,1.3)+int(-1.42*x+3.7,1.3,2.2))
xa5=(int(1.4*x^2,0,0.35)+int(-1.4*x^2+0.98*x,0.35,0.7))/(int(1.4*x,0,
0.35)+int(-1.4*x+0.98,0.35,0.7))
ya5=(int(((x-0.98)/-1.4-x/1.4)*x,0,0.49))/(int(1.4*x,0,0.35)+int(-1.4
*x+0.98,0.35,0.7))
xa6=(int(7.3*x^2,0,0.3)+int(-1.23*x^2+2.56*x,0.3,1.1)+int(-2.43*x^2+3
.88*x,1.1,1.6))/(int(7.3*x,0,0.3)+int(-1.23*x+2.56,0.3,1.1)+int(-2.43
*x+3.88,1.1,1.6))
ya6=(int(((x+1.32)/1.2-x/7.3)*x,0,0.6)+int(((x-2.56)/-1.23-x/7.3)*x,0
.6,2.2))/(int(7.3*x,0,0.3)+int(-1.23*x+2.56,0.3,1.1)+int(-2.43*x+3.88
,1.1,1.6))
xa8=(int(1.24*x^2,0,0.8)+int(0.24*x^2+0.8*x,0.8,5.5)+int(-7.06*x^2+41
.1*x,5.5,5.8))/(int(1.24*x,0,0.8)+int(0.24*x+0.8,0.8,5.5)+int(-7.06*x
+41.1,5.5,5.8))
ya8=(int(((x+40.3)/7.3-(x-0.8)/-1)*x,0,0.8)+int(((x+40.3)/7.3-(x-0.8)
/0.24)*x,0.8,2.2))/(int(1.24*x,0,0.8)+int(0.24*x+0.8,0.8,5.5)+int(-7.
06*x+41.1,5.5,5.8))
xb1=(int(2.19*x^2,0,0.6)+int(1.36*x^2+0.5*x,0.6,4.7)+int(-0.22*x^2+7.93*x,4.7,11.9)+int(-3.79*x^2+50.43*x,11.9,13.3))/(int(2.19*x,0,0.6)+i
nt(1.36*x+0.5,0.6,4.7)+int(-0.22*x+7.93,4.7,11.9)+int(-3.79*x+50.43,1
1.9,13.3))
zb1=(int(((x+42.5)/3.57-(x-0.5)/-0.83)*x,0,0.5)+int(((x+42.5)/3.57-(x
-0.5)/1.36)*x,0.5,5)+int(((x-7.93)/-0.22-(x-0.5)/1.36)*x,5,6.9))/(int
(2.19*x,0,0.6)+int(1.36*x+0.5,0.6,4.7)+int(-0.22*x+7.93,4.7,11.9)+int
(-3.79*x+50.43,11.9,13.3))
xb3=(int(1.64*x^2,0,2.3)+int(-0.77*x^2+5.54*x,2.3,7.2))/(int(1.64*x,0
,2.3)+int(-0.77*x+5.54,2.3,7.2))
zb3=(int(((x+3.24)/1.41-(x-2.3)/-1)*x,0,2.3)+int(((x+3.24)/1.41-(x-2.
3)/0.64)*x,2.3,6.9))/(int(1.64*x,0,2.3)+int(-0.77*x+5.54,2.3,7.2))
xb13=(int(1.25*x^2,0,0.4)+int(-0.83*x^2+0.83*x,0.4,1))/(int(1.25*x,0,
0.4)+int(-0.83*x+0.83,0.4,1))
zb13=(int(-2.08*x^2+0.83*x,0,0.5))/(int(1.25*x,0,0.4)+int(-0.83*x+0.8
3,0.4,1))
xb5=(int(1.29*x^2,0,1.2)+int(-0.88*x^2+2.6*x,1.2,2)+int(-2.13*x^2+5.1
*x,2,2.4))/(int(1.29*x,0,1.2)+int(-0.88*x+2.6,1.2,2)+int(-2.13*x+5.1,
2,2.4))
zb5=(int(((x+2.5)/1.25-(x-2.6)/-2.17)*x,0,0.5)+int((x-2.6)/-0.88-(x-2
.6)/-2.17,0.5,2.6))/(int(1.29*x,0,1.2)+int(-0.88*x+2.6,1.2,2)+int(-2.
13*x+5.1,2,2.4))
xb6=int(-2.17*x^2+2.6*x,0,1.2)/int(-2.17*x+2.6,0,1.2)
zb6=int(((x-2.6)/-2.17)*x,0,2.6)/int(-2.17*x+2.6,0,1.2)
xb8=int(0.55*x^2-0.44*x,0.8,5.5)/int(0.55*x-0.44,0.8,5.5)
zb8=int(((x+0.44)/0.55)*x,0,2.6)/int(0.55*x-0.44,0.8,5.5)
%（5.5.21）式子的解析
j=linspace(0,pi/2,500);
m=linspace(0.004366,0.004366,500);
g=linspace(9.8,9.8,500);
v0=linspace(25,25,500);
H0=linspace(1.75,1.75,500);
k=linspace(0.0007,0.0007,500);
t2=linspace(0,4,500);
t1=m./k.*log(1+k.*v0.*sin(j)./(m.*g));
y1=m./k.*(v0.*sin(j)+m.*g./k).*(1-exp(-k.*t1./m))-m.*g.*t1./k+H0
y2=m.^2.*g./k.^2.*exp(-k.*t2./m)+m.*g.*t2./k-m.^2.*g./k.^2
[AX,H1,H2]=plotyy(y1,j,y2,t2);
set(AX(1),'yTick',[0:0.1:1.6]);
set(get(AX(1),'ylabel'),'string','仰角{theta}(弧度制)');
set(get(AX(2),'ylabel'),'string','下落时间\s')
set(AX(2),'yTick',[0:0.5:4]);
set(AX,'Xlim',[0,40],'xtick',[0:5:40])
legend('等式右边','等式坐标',2)
title('对(5.2.21) 式子的数值解析;k=0.00017')
xlabel('下落高度\m')
%抛投仰角与距离的的解析
j=linspace(0,pi/2,500);
x=linspace(0,pi/2,500);
m=linspace(0.004366,0.004366,500);
g=linspace(9.8,9.8,500);
v0=linspace(25.0998,25.0998,500);
H0=linspace(1.75,1.75,500);
k=linspace(0.00615,0.00615,500);
t2=linspace(0,5,500);
t1=m./k.*log(1+k.*v0.*sin(j)./(m.*g));
y1=m./k.*(v0.*sin(j)+m.*g./k).*(1-exp(-k.*t1./m))-m.*g.*t1./k+H0;
y2=m.^2.*g./k.^2.*exp(-k.*t2./m)+m.*g.*t2./k-m.^2.*g./k.^2;
aa=zeros(500,500);
y11=meshgrid(y1);
for i=1:500;
aa(:,i)=y2'-y11(:,i);
end
bb=aa.^2;
[C,D]=min(bb);
t22=t2(D);
plot(j,t22)
t=t22;
y=(m.*v0.^2.*sin(x).*cos(x)+9.8*m.*v0.*cos(x).*t)./(9.8.*m+k.*v0.*sin
(x));
a=max(y);
b=find(y==a);
plot(x,y);
text(x(b),y(b),['(',num2str(x(b)),',',num2str(y(b)),')'],'color','r')
title('不同k值下的L（{\theta}）函数及最大值点;s=15m');
xlabel('仰角 {\theta}');
ylabel('Suzanne水平运动距离/m');

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可