This example shows how to model a bouncing ball, which is a classical hybrid dynamic system. This model includes both continuous dynamics and discrete transitions. It uses the Symbolic Math Toolbox to help explain some of the theory behind ODE solving in the Simulink® Model of a Bouncing Ball.

The ball rebounds with no angle

There is no drag

Height at time t=0 is 10 m

Thrown with upwards velocity of 15 m/s

The coefficient of restitution is defined as

$${C}_{r}={v}_{b}-{v}_{a}/{u}_{a}-{u}_{b}$$

where v is the speed of object before impact and u is the speed after.

We split the second order differential equation

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

into

$$\frac{dh}{dt}=v$$ discretized as $$\frac{h(t+\Delta t)-h(t)}{\Delta t}=v(t)$$

and

$$\frac{dv}{dt}=-g$$ discretized as $$\frac{v(t+\Delta t)-v(t)}{\Delta t}=-g$$

We will use basic 1st order numerical integration using forward Euler.

$$h(t+\Delta t)=h(t)+v(t)\Delta t$$

$$v(t+\Delta t)=v(t)-g\Delta t$$

Using the Symbolic Math Toolbox, we can approach the problem analytically. This allows us to solve certain problems more easily, such as determining precisely when the ball first hits the ground (below).

Declare our symbolic variables.

syms g t H(t) h_0 v_0

Split the second order differential equation $$\frac{{d}^{2}h}{d{t}^{2}}=-g$$into $$\frac{dh}{dt}=v$$ and $$\frac{dv}{dt}=-g$$.

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

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

Solve the ODE using `dsolve`

:

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

eqn =$$\frac{g\hspace{0.17em}{t}^{2}}{2}+{v}_{0}\hspace{0.17em}t+{h}_{0}$$

Parametrically explore the parabolic profile of motion using `subs`

:

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

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

Find the time at which the ball hits the ground by solving for zero:

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

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

Visualize the solution

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

Format your exact results using variable precision arithmetic with `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 Simulation Parameters

Properties of the ball

`c_bounce = .9; % Bouncing's coefficient of restitution; perfect restitution would be 1`

Properties of the simulation

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)

Declaring the simulation timestep

dt = 0.05; % Animation timestep (s) t_final = 25; % Simulate period (s) t = 0:dt:t_final; % Timespan N = length(t); % Number of iterations

Initialize simulation quantities

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;

Simulate the bouncing ball (we will use basic 1st order numerical integration using forward Euler):

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

Visualize and validate the simulation

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')
```

As a reminder the time of impact was:

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.

From the numerical simulation we can find the closest value in the simulation when $\mathit{h}\left({\mathit{t}}_{\mathit{i}}\right)\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