Этот пример показывает, как смоделировать прыгающий мяч, который является классической гибридной динамической системой. Эта модель включает как непрерывную динамику, так и дискретные переходы. Он использует Symbolic Math Toolbox, чтобы помочь объяснить некоторые теории, лежащие в основе решения ОДУ в Simulink ® Model of a Bouncing Ball .
Мяч отскакивает без угла
Перетаскивание отсутствует
Высота в момент t = 0 составляет 10 м
Отбрасывается со скоростью вверх 15 м/с
Коэффициент реституции определяется как
где v - скорость объекта перед влиянием, а u - скорость после.
Разделим дифференциальное уравнение второго порядка
в
дискретизируется как
и
дискретизируется как
Мы будем использовать базовое численное интегрирование 1-го порядка с помощью forward Euler.
Используя Symbolic Math Toolbox, мы можем подойти к задаче аналитически. Это позволяет нам легче решать определенные задачи, такие как определение именно того, когда мяч впервые ударяется о землю (ниже).
Объявите наши символические переменные.
syms g t H(t) h_0 v_0
Разделите дифференциальное уравнение второго порядка в и .
Dh = diff(H); D2h = diff(H, 2) == g
D2h(t) =
Решить ОДУ используя dsolve
:
eqn = dsolve(D2h, H(0) == h_0, Dh(0) == v_0)
eqn =
Параметрически исследуйте параболический профиль движения с помощью subs
:
eqn = subs(eqn, [h_0, v_0, g], [10, 15, -9.81])
eqn =
Найдите время, в которое мяч попадает в землю, решив на нуль:
assume(t > 0) tHit = solve(eqn == 0)
tHit =
Визуализация решения
fplot(eqn,[0 10]) ylim([0 25])
Форматируйте свои точные результаты с помощью арифметики переменной точности с vpa
:
disp(['The ball with an initial height of 10m and velocity of 15m/s will hit the ground at ' char(vpa(tHit, 4)) ' seconds.'])
The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 3.621 seconds.
Setup параметров симуляции
Свойства мяча
c_bounce = .9; % Bouncing's coefficient of restitution; perfect restitution would be 1
Свойства симуляции
gravity = 9.8; % Gravity's acceleration (m/s) height_0 = 10; % Initial height at time t=0 (m) velocity_0=15; % Initial velocity at time t=0 (m/s)
Объявление временного интервала симуляции
dt = 0.05; % Animation timestep (s) t_final = 25; % Simulate period (s) t = 0:dt:t_final; % Timespan N = length(t); % Number of iterations
Инициализируйте величины симуляции
h=[]; % Height of the ball as a function of time (m) v=[]; % Velocity of the ball (m/sec) as a function of time (m/s) h(1)=height_0; v(1)=velocity_0;
Симулируйте прыгающий мяч (мы будем использовать базовое численное интегрирование 1-го порядка с помощью прямого Эйлера):
for i=1:N-1 v(i+1)=v(i)-gravity*dt; h(i+1)=h(i)+v(i)*dt; % When the ball bounces (the height is less than 0), % reverse the velocity and recalculate the position. % Using the coefficient of restitution if h(i+1) < 0 v(i)=-v(i)*c_bounce; v(i+1)=v(i)-gravity*dt; h(i+1)=0+v(i)*dt; end end
Визуализация и валидация симуляции
plot(t,h,'o') hold on fplot(eqn,[0 10]) title('Height over time') ylim([0 25]) hold off
plot(t,v)
title('Velocity over time')
Напоминаем, что время влияния было:
disp(['The ball with an initial height of 10m and velocity of 15m/s will hit the ground at ' char(vpa(tHit, 4)) ' seconds.'])
The ball with an initial height of 10m and velocity of 15m/s will hit the ground at 3.621 seconds.
Из численного моделирования мы можем найти самое близкое значение в симуляции, когда
i = ceil(double(tHit/dt)); t([i-1 i i+1])
ans = 1×3
3.5500 3.6000 3.6500
plot(t,h,'o') hold on fplot(eqn,[0 10]) plot(t([i-1 i i+1 ]),h([i-1 i i+1 ]),'*R') title('Height over time') xlim([0 5]) ylim([0 25]) hold off