在MATLAB中实现模拟小球上抛和反弹运动

ulppknot

在高中物理中我就学到物体的平抛运动，今天我们在这里也老调重弹下，再次来回顾下这个经典的过程，不过这次讲解的要比之前平抛运动稍微复杂一些：
0 Y( r7 z! ^$ y7 n
/ d6 V0 z( r6 o小球的上抛运动的完整过程描述如下：
2 t0 R2 Q$ F/ N3 h8 ~% b* _
1、我们站在高度为H的看台上抛一个小球
' J! @% t, b5 T3 X6 C) N/ k! O0 [2、上抛的初速度为V，角度为θ
3、小球与空气摩擦力与速度成正比，摩擦系数为μ
4、小球撞地之后能量损失，速度变为原来的k倍，然后继续反弹
( g2 n* L5 W9 W
, g! U+ q3 ?3 w; ^在上面基础上我们继续讨论，比如我们高炮部队，要给予敌方阵地毁灭一击，那么：
/ {! T6 }+ g9 g0 J, Z+ o
1、在已知炮弹初速度V的情况下，以什么角度θ，能使炮弹飞的最远
5 i# d$ s" x# n1 @2、在已知炮弹初速度V和敌方阵地水平距离L的情况下，发射角度θ等于多少时，能整好命中敌方

9 J: `* t( }0 p8 O2 h哈哈，这几个问题好像不是高中的平抛运动能够解决的哦。其实原理和小球的运动方程很容易建立，但是求解起来是有些麻烦的。
& J4 V5 q6 O7 r* e! w8 ]  y/ w9 L* b
我们本次教程这里介绍如何使用MATLAB求解并模拟这个问题，主要设计的内容有：
" V8 K' m" `6 v. O
" X* Y# ~' B  S1 |; U1、小球运动微分法方程求解
0 s8 ]5 e; g$ j. v2、小球着陆时，过零点检测（重点）
7 Z  j; y: p- |% E3 B: d3、小球反弹运动轨迹模拟
, g. @/ C4 P: d! ~7 @  i4、炮弹飞行距离目标最优化
5、发射角度θ的数值求解

4 C7 c2 l% Z* k6 ]: s% j) L* z& Y  A
% W0 m# i3 x9 B由于教程的内容很多，这次主要讲解前三个问题，后两个问题留在稍候的下次教程中讲解！！！


- J* S2 C% y# |% ~9 M# G3 h
小球空中运动方程，只要稍微有一点高中物理和高等数学基础的朋友应该都可以看懂这个方程吧：


/ W: ~3 k, S2 f
上面是相当简单的一个常微分方程组，要求解这个方程组也相当容易，MATLAB提供的ode45函数足以胜任。

, c# S) W# X/ D( n编写主函数main.m

• function main()
• global m mu g
• m=1; % 小球质量
• mu=0.3; % 摩擦系数
• V0=100; % 初始速度
• theta=pi/6; % 初始角度
• g=9.8; % 重力加速
• k=0.8; % 撞击系数
• H=10; % 初始高度
• tspan=[0 10]; % 求解时间范围
• x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值
• [t,x]=ode45(@odefun,tspan,x0);
• X=x(:,1);
• Vx=x(:,2);
• Y=x(:,3);
• Vy=x(:,4);
• figure
• plot(t,x)
• legend('X','Vx','Y','Vy')
• xlabel('时间')
• ylabel('变量')
• grid on
• figure
• plot(X,Y)
• xlabel('水平方向')
• ylabel('竖直方向')
• title('小球运动轨迹')
• grid on
  + T) z! g* |  F' f% f

继续编写微分方程组描述函数odefun.m
% I2 G  M7 u' z
: g7 r' A% H4 D! ]/ D/ \  E

• function dx=odefun(t,x)
• % 微分方程描述函数
• global m mu g
• % X=x(1);
• Vx=x(2);
• % Y=x(3);
• Vy=x(4);
• % 运动微分方程
• dx=[Vx
•     -mu*Vx/m       % dVx/dt
•     Vy
• -mu*Vy/m-g];   % dVy/dt
  

7 A, G$ g( s6 }. ~5 h+ G8 M1 H0 z; C
0 z! U. Q/ m9 p8 E/ W- q' y7 S运行之后我们得到如下结果：

. |" g5 r  l8 Q/ N) k


5 V8 L9 D0 M+ j( _8 t从上面的“小球运动轨迹图形”可看出，小球在水平位移265米左右就撞地了，小球撞地以后会反弹，显然上图撞地之后的数据就是无效的。

但是要模拟小球的着陆撞击反弹过程，我们希望程序能够自动检测到那个小球撞击的时刻，听起来好像有些深奥哦，其实这个也没什么高深的，因为MATLAB的ode45函数早就为我们考虑到这点了，此时只需要设置ODE方程的一个过零点检测，也就是设置ode45函数的events属性即可。

修改main.m函数，主要是通过options参数添加events事件检测（黄色部分）：
2 ~! t& p2 t3 v

• function main
• global m mu g
• m=1; % 小球质量
• mu=0.3; % 摩擦系数
• V0=100; % 初始速度
• theta=pi/6; % 初始角度
• g=9.8; % 重力加速
• k=0.8; % 撞击系数
• H=10; % 初始高度
• tspan=[0 10]; % 求解时间范围
• x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值
• options=odeset('events',@events); % 设置小球过零点（撞地）检测条件
• [t,x,te,xe,ie]=ode45(@odefun,tspan,x0,options);
• X=x(:,1);
• Vx=x(:,2);
• Y=x(:,3);
• Vy=x(:,4);
• figure
• plot(t,x)
• legend('X','Vx','Y','Vy')
• xlabel('时间')
• ylabel('变量')
• title('参数随时间变化图形')
• grid on
• figure
• plot(X,Y)
• xlabel('水平方向')
• ylabel('竖直方向')
• title('小球运动轨迹图形')
• grid on
  ( I; x, R  X, W7 S; M

' N; P+ O- K3 G7 w$ \8 p

编写events.m文件描述需要检测的事件（过零点检测）：
3 |+ Z- B3 c" @+ r( B& v  A

• function [value,isterminal,direction]=events(t,x)
• % 事件检查函数，此时需要做的是过零点检测
• Y=x(3);
• % ode45函数自动检查当value=0是否成立
• % 在此问题上我们要求检测Y=0的点，因此设置value=Y
• % 如果我们要检测Y=2，那么就设置value=Y-2
• value=Y;
• % 检测到指定条件时，是否终止ode45函数的运行
• % 1表示终止，0表示继续
• % 在我们这个问题上，我们只要检测到零点时就停止程序
• isterminal=1;
• % value过零点检测的方向
• % -1表示由正到负，+1表示由负到正
• % 对于我们这个问题，当然是由正到负
• direction=-1;
  

. I! f" x# Q8 n8 r
5 d3 t2 B/ O  j3 I& u+ G* X
再次运行程序，得到如下运动轨迹图：
/ ?3 h7 l0 R( g4 p
, o  ?( P5 S/ x( M" ~8 ?

从图形可以看出，程序在检测到小球着陆，也就是Y=0时自动终止程序运行了，并通过ode45的te和xe参数返回事件发生的时刻te和此时所有的状态参数xe。

之前的工作我们已经完成了方程求解和零点检测，但是还是没有完整的描述整个小球上抛运动的过程呀，比如反弹就没有。我们希望程序能够自动将整个过程都能够展现出来。
( f. N9 j2 ]8 |9 R% g- ?+ ]0 k
有些网络是不是想说：“是不是MATLAB的ode45函数也为我们考虑到这点，并早已经准备好了某个函数或者参数提供给我们直接调用呀？”哈哈，这次很惋惜的告诉您，这个需要我们自己动手丰衣足食了，MATLAB不总是靠得住的，毕竟软件是死的，人是活的。。。。

, d) l. v( `4 \) e- F4 O对于反弹以及全过程描述，其实也不难！我们只要在检测到撞地点以后，以撞击时候的参数为初值继续求解微分方程，这样一直下去。最后将所有的小过程连接起来，于是完整的小球上抛和反弹过程就可以展现出来了！每个小过程的求解，我们可以在循环语句中进行！

' U% j! `  `+ j8 l4 Q& ^4 e) C于是重新修改（黄色部分）main.m函数：

• function main
• global m mu g
• m=1; % 小球质量
• mu=0.25; % 摩擦系数
• V0=100; % 初始速度
• theta=pi/6; % 初始角度
• g=9.8; % 重力加速
• k=0.9; % 撞击系数
• H=10; % 初始高度
• tspan=[0 100]; % 求解时间范围
• x0=[0,V0*cos(theta),H,V0*sin(theta)]; % 微分方程初值
• options=odeset('events',@events); % 设置小球过零点（撞地）检测条件
• n=10; % 反弹次数
• % 初始化必要的参数
• tout=tspan(1); % 记录完整过程的时间
• xout=x0; % 记录完整过程的参数
• teout=[]; % 记录撞击点时刻
• xeout=[]; % 记录撞击点参数
• 
    |% V  h' {5 C, h
• for i=1:n
•     [t,x,te,xe]=ode45(@odefun,tspan,x0,options);
•     % 保存该阶段的参数
•     tout=[tout;t(2:end)];
•     xout=[xout;x(2:end,: )];
•     teout=[teout;te];
•     xeout=[xeout;xe];
•    
•     % 修改方程求解区间
•     tspan(1)=te;
•     % 修改初值方程
•     x0(1)=xe(1); % 水平位移
•     x0(2)=xe(2)*k; % 水平速度，注意乘以撞击系数
•     x0(3)=0; % 竖直位移
•     x0(4)=-xe(4)*k; % 竖直速度，注意乘以撞击系数，并且反弹
• end
• X=xout(:,1);
• Vx=xout(:,2);
• Y=xout(:,3);
• Vy=xout(:,4);
• figure
• plot(tout,xout)
• hold on
• plot(teout,xeout,'ro','MarkeRFaceColor','red')
• legend('X','Vx','Y','Vy')
• xlabel('时间')
• ylabel('变量')
• title('参数随时间变化图形')
• grid on
• figure
• plot(X,Y)
• hold on
• plot(xeout(:,1),xeout(:,3),'ro','MarkerFaceColor','red')
• xlabel('水平方向')
• ylabel('竖直方向')
• title('小球运动轨迹图形')
• grid on
  ; h8 X6 C9 `7 W$ A/ k

0 q7 d  b. G3 `- V
2 K' f  Z; ~+ H# p
重新运行main函数，可以到小球上抛以反弹运动的完整过程了：


17 分钟前 上传
! W) c' x/ {7 K& t# n# a; Y: g下载附件 (36.32 KB)
1 t- C/ N6 _( H) `# S

4 Z" U  D. w- y6 P7 m' w
7 g2 ]( p% q; i
) I  F1 v7 N* N8 {8 M2 n

CCxiaom

在MATLAB中实现模拟小球上抛和反弹运动