Проверьте модель Simulink с помощью Symbolic Math Toolbox

Этот пример показывает, как смоделировать прыгающий мяч, который является классической гибридной динамической системой. Эта модель включает как непрерывную динамику, так и дискретные переходы. Он использует Symbolic Math Toolbox, чтобы помочь объяснить некоторые теории, лежащие в основе решения ОДУ в Simulink ® Model of a Bouncing Ball .

Предположения

  • Мяч отскакивает без угла

  • Перетаскивание отсутствует

  • Высота в момент t = 0 составляет 10 м

  • Отбрасывается со скоростью вверх 15 м/с

Происхождение

Коэффициент реституции определяется как

Cr=vb-va/ua-ub

где v - скорость объекта перед влиянием, а u - скорость после.

Разделим дифференциальное уравнение второго порядка

d2hdt2=-g

в

dhdt=v дискретизируется как h(t+Δt)-h(t)Δt=v(t)

и

dvdt=-g дискретизируется как v(t+Δt)-v(t)Δt=-g

Мы будем использовать базовое численное интегрирование 1-го порядка с помощью forward Euler.

h(t+Δt)=h(t)+v(t)Δt

v(t+Δt)=v(t)-gΔt

Решите задачу аналитически

Используя Symbolic Math Toolbox, мы можем подойти к задаче аналитически. Это позволяет нам легче решать определенные задачи, такие как определение именно того, когда мяч впервые ударяется о землю (ниже).

Объявите наши символические переменные.

syms g t H(t) h_0 v_0

Разделите дифференциальное уравнение второго порядка d2hdt2=-gв dhdt=v и dvdt=-g.

Dh = diff(H);
D2h = diff(H, 2) == g
D2h(t) = 

2t2 H(t)=gdiff (H (t), t, 2) = = g

Решить ОДУ используя dsolve:

eqn = dsolve(D2h, H(0) == h_0, Dh(0) == v_0)
eqn = 

gt22+v0t+h0(g * t ^ 2 )/2 + v_0*t + h_0

Параметрически исследуйте параболический профиль движения с помощью subs:

eqn = subs(eqn, [h_0, v_0, g], [10, 15, -9.81])
eqn = 

-981t2200+15t+10- (981 * t ^ 2 )/200 + 15 * t + 10

Найдите время, в которое мяч попадает в землю, решив на нуль:

assume(t > 0)
tHit = solve(eqn == 0)
tHit = 

20526109+500327(20 * sqrt (sym (5)) * sqrt (sym (26)) )/109 + sym (500/327)

Визуализация решения

fplot(eqn,[0 10])
ylim([0 25])

Figure contains an axes. The axes contains an object of type functionline.

Форматируйте свои точные результаты с помощью арифметики переменной точности с 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

Figure contains an axes. The axes with title Height over time contains 2 objects of type line, functionline.

plot(t,v)
title('Velocity over time')

Figure contains an axes. The axes with title Velocity over time contains an object of type line.

Валидация цифр с помощью аналитики

Сравните свои аналитические результаты с числовыми результатами.

Напоминаем, что время влияния было:

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.

Из численного моделирования мы можем найти самое близкое значение в симуляции, когда h(ti)0

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

Figure contains an axes. The axes with title Height over time contains 3 objects of type line, functionline.