3,Matlab仿真弹跳球

原理

使用ode45的event函数，对运动方程积分，遇到特定实践（碰撞）则进行改变初值条件，重新开始积分

matlab代码

运动方程

小球在空中的运动方程为

function dqdt = rhs(t,y,ball)
% y = [y ydot]
ydot = y(2);
yddot = -ball.g;

dqdt = [ydot yddot]';

onebounce函数

可以将小球的弹跳分为多次过程，用循环实现，每次弹跳的代码基本相同

function [tt,yy] = my_bounce(t0,ystart,ball)

dt = 10; % 随便给，需要大一点
%%%% 计算一次碰撞前的t和y向量
options = odeset('AbsTol',1e-8,'RelTol',1e-8,'Events',@mycontact);
[t y] = ode45(@rhs,[t0 t0+dt],ystart,options,ball); %检测Event发生则积分终止

% ynew(1,1) = y(end,1);
y(end,2) = -ball.e*y(end,2);%终止后，小球回弹，设置反弹的速度

tt = t;
yy = y;

做动画函数

为了直观显示，通常将仿真结果做成动画，代码如下：

function myanimation(t_all,y_all,ball)
tstart = t_all(1);
tend = t_all(end);
tinter = linspace(tstart,tend,ball.movieFps*(tend-tstart));% 建立每一帧的时间节点

[m,n] = size(y_all);

for i = 1:n
    yinter(:,i) = interp1(t_all,y_all(:,i),tinter);% 使用插值法，将每个时间节点都对应上相应的数据
end

for i = 1:length(tinter)
    plot(0,yinter(i,1),'ro','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',20); hold on;
    line([-0.1 0.1],[ball.ground_animate ball.ground_animate],'LineWidth',2);
    axis([-0.1 0.1 min(yinter(:,1)-1) max(yinter(:,1)+1)]);
    pause(ball.moviePause);
    hold off;
end

mycontact函数

function [gstop,isterminal,direction] = mycontact(t,y,ball)

gstop = y(1) - ball.ground;
direction = -1;%往下走的过程中，检测为碰撞
isterminal = 1;

main部分

clc
clear
close all
%%%% 建立对象球的参数 %%%%
ball.g = 9.8;
ball.e = 0.9;%反弹系数
ball.ground = 0;
ball.ground_animate = -0.2;%做动画用的地面

ball.movieFps = 15;
ball.moviePause = 0.05;

%%%% 初始条件 %%%%
y0 = 2; y0dot = 10;
ystart = [y0 y0dot];

t0 = 0; tend = 10;
tstart = t0;
%%%% 仿真计算 %%%%
t_all = tstart; y_all = ystart;
while(t0<tend)
    [t y] = my_bounce(t0,ystart,ball); %仿真一次计算
    t_all = [t_all;t(2:end)];% 把时间向量记录进来,从2开始的原因是t(1)就是上次循环的t(end)
    y_all = [y_all;y(2:end,:)];

    t0 = t(end);
    ystart = y(end,:);
end

figure(1)
subplot(2,1,1);
plot(t_all,y_all(:,1),'r','LineWidth',2);
ylabel('纵方向位置','FontSize',12);
xlabel('时间','FontSize',12);
subplot(2,1,2);
plot(t_all,y_all(:,2),'b','LineWidth',2);
ylabel('纵方向速度','FontSize',12);
xlabel('时间','FontSize',12);

figure(2)
myanimation(t_all,y_all,ball);

仿真结果