Подтвердите модель Simulink Используя Symbolic Math Toolbox

Скрипт Open Live Script

В этом примере показано, как смоделировать прыгающий мяч, который является классической гибридной динамической системой. Эта модель включает и непрерывную динамику и дискретные переходы. Это использует Symbolic Math Toolbox, чтобы помочь объяснить часть теории позади решения ОДУ в Модели Simulink® Прыгающего мяча.

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

Мяч оживляется без угла
Нет никаких, перетаскивают
Высота во время t=0 составляет 10 м
Выданный с вверх скоростью 15 м/с

Деривация

Коэффициент восстановления задан как

$C_{r} = v_{b} - v_{a} / u_{a} - u_{b}$

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

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

$\frac{d^{2} h}{d t^{2}} = - g$

$\frac{d h}{d t} = v$ дискретизированный как $\frac{h (t + Δ t) - h (t)}{Δ t} = v (t)$

$\frac{d v}{d t} = - g$ дискретизированный как $\frac{v (t + Δ t) - v (t)}{Δ t} = - g$

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

$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

Разделите дифференциальное уравнение второго порядка $\frac{d^{2} h}{d t^{2}} = - g$ в $\frac{d h}{d t} = v$ и $\frac{d v}{d t} = - g$ .

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

D2h(t) = 
 $\frac{\partial^{2}}{\partial t^{2}} H (t) = g$

Решите ОДУ с помощью dsolve:

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

eqn = 
 $\frac{g t^{2}}{2} + v_{0} t + h_{0}$

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

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

eqn = 
 $- \frac{981 t^{2}}{200} + 15 t + 10$

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

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

tHit = 
 $\frac{20 \sqrt{5} \sqrt{26}}{109} + \frac{500}{327}$

Визуализируйте решение

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

Figure contains an axes object. The axes object 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 object. The axes object with title Height over time contains 2 objects of type line, functionline.

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

Figure contains an axes object. The axes object 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 (t_{i}) \approx 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 object. The axes object with title Height over time contains 3 objects of type line, functionline.

Документация