MATLAB tutorial for the first course, part 1.2: Applications

========================== to be checked ========================

Vertical Motion

It is known from elementary physics that, in the absence of air friction, a ball thrown up from the ground into earth's atmosphere with initial speed v₀ would attain a maximum altitude of $ v_0^2 /(2g) .$ In this case the return time is $ 2\, v_0 /g, $ independent of the ball's mass. Here g is the acceleration due to gravity. If the ball is thrown up from altitude y₀ (which we later assume to be zero), then the time T₀ spent traveling is given by

\[ \left\{ \begin{array}{c} \mbox{travel time with no air resistance} \\ \mbox{when thrown from height $y_0$} \end{array} \right\} \ = \ T_0 \ =\ \frac{v_0 + \sqrt{v_0^2 +2\, y_0\, g}}{g} \]

The presence of air influences the ball's motion: it experiences two forces acting on it---the force of gravity and the air resistance force. Let's define the symbol T as follows:

\[ \left\{ \begin{array}{c} \mbox{travel time with air resistance} \\ \mbox{when thrown from height $y_0$} \end{array} \right\} \ = \ T \]

Without air resistance, the object travels farther up than with air resistance. On the way down, without air resistance the object travels a larger distance, but it also gathers more speed. A natural question is, which travel time (with air resistance vs. no air resistance) is larger? Also, it is of interest to find the maximum altitude $ y_{\max} $ of the ball, the time T_max to reach maximum altitude, and the time $ T_{\text{down}} $ to return back from y_max. Therefore, $ T_{\max} + T_{\text{down}} = T_{\text{total}} $ (the total time the ball spent in the air). The landing velocity is denoted by $ v_{\ell} . $

Example 1: Air resistance is the force that acts in the direction opposite to the motion of an object through air. Air resistance depends on the shape, material, and orientation of the object, the density of the air, and the object's relative speed.

We would like to think that there is a nice formula for the air resistance in terms of speed and other variables. Such a formula would help in making calculations and predicting various quantities. A starting point for obtaining such a formula is our everyday experience. Based on our experience, a reasonable assumption to make\footnote{Our intuition based on everyday experience is limited to a small range of conditions. This may lead to erroneous assumptions.

It has been observed that, under suitable conditions, the magnitude of the air resistance is proportional to a power of the speed s=|v|:

\begin{equation} \tag{M} F \propto s^p, \qquad \mbox{which may be written as}\qquad F(v) = k \, |v|^p. \end{equation}

Here v is velocity, and both k and p are positive constants. For very small objects, such as a speck of dust (about 1 micrometer or 0.001mm), p=1 seems to give a reasonable formula for the air resistance. For larger, human scale objects moving at relatively large speed, p=2 works better. Therefore, the magnitude of the air resistance F as a function of velocity v is assumed be given by formula (M).

The air resistance force depends on the velocity (v) of the object at time t, so let us denote this force with the symbol F(v). Note that the air resistance, force F(v), always acts in the direction opposite to the motion. Therefore, F(v) acts in the down (negative) direction when the ball is moving up, and it acts in the up (positive) direction when the ball is moving down. If we measure the displacement y = y(t) vertically upwards from the ground, then $ v= {\text d}y/{\text d}t = \dot{y} $ is the velocity of the object. Newton's law of motion for the ball on the way up gives the differential equation

\begin{equation} \label{vert.1} m\, \dot{v} = -m\,g -F(v) \ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g - \frac{1}{m}\, F(v), \end{equation}

and on the way down,

\begin{equation} \label{vert.2} m\, \dot{v} = -m\,g +F(v)\ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g + \frac{1}{m}\, F(v). \end{equation}

Since we assume $ F(v) = k \, |v|^p , $ the equation of motion on way up becomes

\begin{equation} \label{vert.3} m\, \dot{v} = -m\,g -k\,|v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p , \end{equation}

and on the way down,

\begin{equation} \label{vert.4} m\, \dot{v} = -m\, g +k\, |v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p . \end{equation}

To find an equation for $ v_{\ell} , $ the landing velocity, we rewrite the equation $ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p $ as $ {\text d}\,t = - {\text d}v/(g +F(v)/m) $ and integrate both sides from t=0 and v = v₀ to $ t=T_{\max} $ and v=0. Here $ T_{\max} $ is the time to reach the maximum altitude $ y_{max} , $ which is also the time to have velocity v=0. We obtain,

\begin{equation} \label{Tmax.1} \int_{0}^{T_{max}}\,{\text d}t = - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} \, , \quad \mbox{that is,} \quad T_{max} = - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} \end{equation}

A similar formula is valid for time T. From $ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p , $ we get

\begin{equation} \label{Tmax.2} \int_{T_{max}}^T \, {\text d}t = \int_0^{v_\ell} \frac{ \,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \quad\Longrightarrow \quad T - T_{\max} = \int_0^{v_\ell} \frac{ \,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \end{equation}

Equating the time $ T_{max} $ in the above equations, we have

\begin{equation} \label{Tlanding} T \ = \ - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} + \int_0^{v_\ell} \frac{{\text d}\,v}{-g+\frac{k}{m}\,|v|^p} \end{equation}

This equation gives T as a function of $ v_\ell . $ The next step is to find an equation for $ v_\ell . $

To find an equation for $ v_\ell , $ we rewrite the equation $ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p $ as $ v\,{\text d}t = -v\, {\text d}v/(g +F(v)/m) $ and integrate both sides from t=0 and v = v₀ to $ t=T_{\max} $ and v=0. Here $ T_{\max} $ is the time to reach the maximum altitude $ y_{\max} , $ which is also the time to have velocity v=0. Using the fact that the integral of the velocity is the displacement, we obtain,

\begin{equation} \label{ymax.1} \int_{0}^{T_{max}} v\,{\text d}t = - \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \, , \quad \mbox{that is,} \quad y_{max} - y_0 = - \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \end{equation}

A similar formula is valid for the distance traveled down. From equation $ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p , $ we get

\begin{equation} \label{ymax.2} \int_{T_{max}}^T v\, {\text d}t = \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \quad\Longrightarrow \quad - y_{\max} = \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \end{equation}

Equating the maximum distance $ y_{\max} $ traveled in both ways, we get an equation involving the landing velocity $ v_{\ell} : $

\begin{equation} \label{vlanding} \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \ = \ - y_0 + \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \end{equation}

This equation is an equation where the unknown is $ v_{\ell} , $ which does not appear explicitly solved for.

g := 9.806
Solve[Integrate[v/(-g + 0.01*v^2), {v, 0, vl}] == Integrate[v/(g + 0.01*v^2), {v, 50, 0}], vl]

Out[2]= {{vl -> -26.5393}, {vl -> 26.5393}}

We would like to determine the ratio:

\[ \gamma = T/T_0 = Tg/ \left( v_0 + \sqrt{v_0^2 +2 x_0 g}\right), \]

where T time in air with air resistance and T₀ is the time in air without air resistance.

For a tennis ball thrown upward with the initial velocity v₀ =10, it is possible to find x₀ that γ > 1 when p=0.9. In general, it is unknown for what values of p< 1, x₀, and v₀ we can achieve γ > 1.

In the code above,

%% Up and Down for a ball with drag and no drag
%% for the simulation, we compare the velocity and trajectory of two balls one with no drag and one with drag, this will be utilized to understand the effects on, velocity, acceleration and height that these two identical balls have with the differing conditions of drag and other factors
%% such an experiment is interesting in the context of physics as it allows us to be able to compare and contrast different effects removed from the
%% factors in nature outside of our control and as such is an interesting %% and worthwile simulation ■

Example 2: We consider a model of falling object, say a tennis ball, to a flat surface that moves up and down periodicaly (from racket of a tennis player). Using vertical axis directed upward, we denote v(t) as the velocity of the ball and y(t) as its position/height at time t. It has been observed that, under suitable conditions, the magnitude of the air resistance is proportional to the power of speeds=|v|:

\begin{equation} \tag{M} F \propto s^p, \qquad \mbox{which may be written as}\qquad F(v) = k \, |v|^p. \end{equation}

Note that the air resistance force F(v) always acts in the direction opposite to the motion. Therefore, F(v) acts in the down (negative) direction when the ball is moving up, and it acts in the up (positive) direction when the ball is moving down.

Suppose that initially at t=0, the ball of mass m is dropped from the altitude $ y=h > 1 $ without initial velocity. At the same time, it is assumed that the floor starts moving according to the formula $ z= \sin \omega t. $ When elastic ball hits a hard flat surface, it bounces back with the same velocity. It is assumed that the collision is totally elastic, so the ball loses no kinetic energy in the collision, and its speed after collision is the same as before the collision. At this point, ignore the time needed for the ball to be deformed during collision before fully rebounded and has lifted off from the surface instantly. Hence the ball can be treated as a rigid body with negligible deformation during impact.

After collision, the ball climbs up until its velocity becomes zero, and then the ball falls vertically downward under the influence of gravity, hits the the moving floor, and bounces back.

Derivation of a differential equation

Newton's law of motion for the ball on the way down is

\begin{equation*} %\label{vert.1} m\, \dot{v} = -m\,g +F(v)\ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g + \frac{1}{m}\, F(v), \end{equation*}

and on the way up

\begin{equation*} %\label{vert.2} m\, \dot{v} = -m\,g -F(v) \ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g - \frac{1}{m}\, F(v), \end{equation*}

where g is the acceleration due to gravity. Since we assume $ F(v) = k \, |v|^p ,$ the equation of motion on the way down becomes

\begin{equation*} %\label{vert.3} m\, \dot{v} = -m\, g +k\, |v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p . \end{equation*}

and on way up

\begin{equation} %\label{vert.4} m\, \dot{v} = -m\,g -k\,|v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p , \end{equation}

Assuming that p=2 for a tennis ball, the above differential equations can be integrated using separation of variables:

\begin{equation} %\label{vert.1} \int \frac{{\text d}v}{-g+\frac{k}{m}\,|v|^2} = \int {\text d}t \qquad\Longrightarrow \qquad \sqrt{\frac{m}{g\,k}}\, \mbox{Arctanh} \left( v\,\sqrt{\frac{k}{g\,m}} \right) = t , \end{equation}

and

\begin{equation} %\label{vert.2} \int \frac{{\text d}v}{g+\frac{k}{m}\,|v|^2} = -\int {\text d}t \qquad\Longrightarrow \qquad \sqrt{\frac{m}{g\,k}}\, \mbox{Arctan} \left( v\,\sqrt{\frac{k}{g\,m}} \right) = -t , \end{equation}

where

\[ \mbox{Arctanh} x = \frac{1}{2} \, \ln \frac{1+x}{1-x} \quad \mbox{for} \quad |x| < 1. \]

Since the velocity v(t) is the derivative of ball's position y(t), you may need to integrate

\[ \int \arctan (\alpha v)\, {\text d}v = v\,\arctan (\alpha v) - \frac{1}{2\alpha} \, \ln \left( 1 + \alpha^2 v^2 \right) , \quad \int \mbox{arctanh} (\alpha v)\, {\text d}v = v\,\mbox{arctanh} (\alpha v) - \frac{1}{2\alpha} \, \ln \left\vert 1 - \alpha^2 v^2 \right\vert . \]

For example, at the initial stage, you need to solve two IVPs:

\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(0) =0; \qquad \frac{{\text d}y}{{\text d}t} = v, \quad y(0) =h . \]

Input parameters

$ g \approx 9.806 $ m/sec²	the acceleration due to gravity near sea level at 45 deg. latitude;
m = 0.08	mass of the object, in kg because a tennis ball is about 80 grams;
k = 0.02	drag coefficient, positive;
p = 2	power of the speed term in the resistance force;
ω = π	frequency of the oscillating floor;
$ y_0 =h > 1 $	initial altitude, positive, in meters.

Derivation of solution

A ball that is dropped from height h> 1 can be described by its velocity v(t) and position y(t):

\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(0) =0; \qquad \frac{{\text d}y}{{\text d}t} = v, \quad y(0) =h . \]

Separation of variables yields

\[ \int_0^v \frac{{\text d}v}{-g+\frac{k}{m}\,|v|^2} = \int_0^t {\text d}t \qquad \Longrightarrow \qquad \mbox{Arctanh} \left( \sqrt{\frac{k}{gm}} \, v \right) = - \sqrt{\frac{gk}{m}} \,t . \]

Therefore, we find the velocity v(t) on first stage of falling from hight h> 1 explicitly:

\[ v(t) = - \sqrt{\frac{gm}{k}} \,\tanh \left( \sqrt{\frac{gk}{m}} \,t \right) , \qquad k\,v^2 < g\,m. \]

Then we find ball's position by integrating v(t):

\[ y(t) = h- \frac{m}{k}\,\ln \cosh \left( \sqrt{\frac{gk}{m}} \,t \right) . \]

This equation is valid untill the ball meets the surface $ z= \sin \pi t. $ Therefore, we need to solve the transcendent equation:

\[ \sin \pi t = h- \frac{m}{k}\,\ln \cosh \left( \sqrt{\frac{gk}{m}} \,t \right) . \]

Mathematica confirms:

Assuming[v > 0 && k > 0 && m > 0 && g > 0, Integrate[1/(-g + k/m*v1^2), {v1, 0, v}]]
Integrate[A*Tanh[B*t], t]
k := 0.02; g := 9.806; m := 0.08; h0 := 2; (* h0 is the initial height *)
FindRoot[Sin[\[Pi] t] == 2 - 4* Log[Cosh[1.5657266683556232*t]], {t, 0.4}]
Plot[{Sin[Pi*t], h[t]}, {t, 0.4, 0.6}]

Since Mathematica provides its approximate value to be 0.4716548296910227, we denote it by T1.

On the second stage, the ball bounced up with the initial velocity $ V1 \approx 3.93453 $ and from the position $ Y1 \approx 0.996038. $ Therefore, we need to solve two intial value problems:

\[ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^2 , \quad v(T1) = V1, \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T1) = Y1, \quad \mbox{for }t\ge T1. \]

Separating variables and integrating, we obtain

\[ v(t) = V1 + \sqrt{\frac{gm}{k}} \, \tan \left[ \sqrt{\frac{gk}{m}} \left( T1 -t \right) \right] , \quad T1 \le t \le T2 , \]

where T2 is the value of time when v(T2) =0. Then we find ball's position:

\[ y(t) = Y1+ V1\left( t- T1 \right) +\frac{m}{k} \,\ln \cos \sqrt{\frac{gk}{m}} \left( T1 -t \right) , \quad T1 \le t \le T2 . \]

Mathematica find its value to be $ T2 \approx 0.829902 : $

v2[t_] := V1 + Sqrt[g*m/k]*Tan[Sqrt[g*k/m]*(T1 - t)]
FindRoot[V1 + v2[t] == 0, {t, 0.8}]
Y2 := y[t] /. t -> T2

The position of the ball at t= T2 will be $ y(T2) \approx 1.74026 . $ After T2, the ball starts falling down and we need to solve the initial value problems:

\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(T2) =0; \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T2) =Y2 \approx 1.74026 . \]

We again separate variables and obtain

\[ v(t) = \sqrt{\frac{gm}{k}} \, \tanh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] , \quad T2 \le t\le T3 , \]

\[ y(t) = Y2 + \frac{m}{k} \, \ln \cosh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] , \quad T2 \le t\le T3 . \]

To find T3, we need to solve the transcendent equation:

\[ \sin \left( \omega t \right) = Y2 + \frac{m}{k} \, \ln \cosh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] , \]

so we ask {\em Mathematica:}

v3[t_] = Sqrt[g*m/k]*Tanh[Sqrt[g*k/m]*(t - T2)]
FindRoot[Sin[\[Pi] t] == Y2 - 4*Log[Cosh[1.5657266683556232*(t - T2)]], {t, 1.5}]
T3 := 1.6458839181239935
y[t_] = Piecewise[{{2 - m/k*Log[Cosh[Sqrt[g*k/m]*t]], 0 < t < T1}, {Y1 + V1*(t - T1) +
m/k*Log[Cos[Sqrt[g*k/m]*(T1 - t)]], T1 < t <= T2}, T1 < t <= T2}, {Y2 - 4*Log[Cosh[Sqrt[g*k/m]*(t - T2)]], T2 < t <= T3}}]
Plot[y[t], {t, 0, T3}]

At $ t= T3 \approx 1.64588 , $ the ball meets the floor and start climbing up. So we its velocity should be the solution of the following initial value problem:

\[ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^2 , \quad v(T3) = V3, \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T3) = Y3, \quad \mbox{for } T3 \le t \le T4, \]

where T4 is time where the velocity is zero, and $ Y3 = y(T3) = \sin \left( \pi T3 \right) . $ Mathematica provides the numerical values:

Y3 = y[t] /. t -> T3
v4[t_] = Sqrt[g*m/k]*Tan[Sqrt[g*k/m]*(T3 - t)]
V3 = v3[t] /. t -> T3
FindRoot[V3 + v4[t] == 0, {t, 2.0}]
T4 := 2.097992138322688

With $ V3 \approx 5.36008 , $ we find velocity

\[ v(t) = V3 + \sqrt{\frac{gm}{k}} \, \tan \left[ \sqrt{\frac{gk}{m}} \left( T3 -t \right) \right] , \quad T3 \le t \le T4 , \]

and position of the ball:

\[ y(t) = Y3+ V3\left( t- T3 \right) +\frac{m}{k} \,\ln \cos \sqrt{\frac{gk}{m}} \left( T3 -t \right) , \quad T3 \le t \le T4 . \]

Starting with t= T4, the ball falls till it meets the floor, which find with the following command:

FindRoot[Sin[\[Pi] t] == Y3 + V3*(t - T3) + m/k*Log[Cos[Sqrt[g*k/m]*(T3 - t)]], {t, 2.1}]
T5 := 2.1360983913469607

Then we solve the initial value problems in the time interval [T5, T6]:

\[ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^2 , \quad v(T5) = V5, \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T5) = Y5, \quad \mbox{for } T5 \le t \le T6, \]

where T6 is time when ball's velocity is zero.

FindRoot[Sin[\[Pi] t] == Y3 + V3*(t - T3) + m/k*Log[Cos[Sqrt[g*k/m]*(T3 - t)]], {t, 2.1}]
T5 := 2.1360983913469607
Y4 := (Y3 + V3*(t - T3) + m/k*Log[Cos[Sqrt[g*k/m]*(T3 - t)]]) /. t -> T4
y[t_] = Piecewise[{{2 - m/k*Log[Cosh[Sqrt[g*k/m]*t]], 0 < t < T1}, {Y1 + V1*(t - T1) +
m/k*Log[Cos[Sqrt[g*k/m]*(T1 - t)]],
T1 < t <= T2}, {Y2 - 4*Log[Cosh[Sqrt[g*k/m]*(t - T2)]],
T2 < t <= T3}, {Y3 + V3*(t - T3) +
m/k*Log[Cos[Sqrt[g*k/m]*(T3 - t)]],
T3 < t <= T4}, {Y4 - 4*Log[Cosh[Sqrt[g*k/m]*(t - T4)]], T4 < t <= T5}}]
Plot[{Sin[Pi*t], y[t]}, {t, 0, T5}]

Graphs of sin π t and position of *y(t)*.

Finally, we use Euler's method to find approximate solution to the sequence of initial value problems.


function jumping
g = 9.807;        % gravitational acceleration, m/s^2
A = 1; w = .5;    % amplitude, m, and frequency of surface oscillation, Hz or 1/s 
m = .08; k = .02; % mass of a ball, kg, and coefficient of air resistance, kg/m
H = 2.5;          % initial height of a ball, m
v0 = 0;           % initial speed of a ball, m/s
J = 5;            % maximal number of bounces
T = 5;            % maximal duration of a process, s
dt = 1e-5;        % time integration step, s 
epsh = .0001;     % a threshold height of a ball over surface, m
epsv = .0005;     % a threshold for speed of a ball, m/s
n = round(T/dt);  % maximal number of steps
s = @(t) A*sin(2*pi*w*t);         % surface oscillation function
vs = @(t) A*2*pi*w*cos(2*pi*w*t); % surface velocity function
i = 1;   % steps counter
j = -1;  % bounces counter
h(n) = 0; v(n) = 0;   % memory allocation
h(1) = H; v(1) = v0;  % initial values

while j <= J
  while i < n % down
    v(i + 1) = (g - k*v(i)^2/m)*dt + v(i);
    h(i + 1) = -v(i)*dt + h(i);
    S = s(dt*(i + 1));    
    if abs(S - h(i + 1)) < epsh % a ball meets a surface
      h(i + 1) = S + epsh;
      v(i + 1) = v(i) + vs(i*dt);
      i = i + 1;
      j = j + 1;
      break
    end   
    i = i + 1;
  end
  
  if j >= J || i >= n
   break 
  end

  while i < n  % up
    v(i + 1) = -(g + k*v(i)^2/m)*dt + v(i);    
    h(i + 1) = v(i)*dt + h(i);
    S = s(dt*(i + 1));  
    if abs(v(i + 1)) < epsv % speed of a ball near zero
      break
    end
     if (abs(S - h(i + 1))) < epsh % a ball meets a surface
       h(i + 1) = S + epsh;
       v(i + 1) = v(i) + vs(i*dt);  
       j = j + 1;
       if j >= J || i >= n  
         break
       end
     end  
    i = i + 1;
  end
end

subplot(2,1,1)
t = dt*(0:i - 1);
plot(t,h(1:i),t,s(dt*(1:i)),'linewidth',2);
title('Height, m');  xlabel('Time, s')
legend('Ball','Surface')
grid; axis tight; 
subplot(2,1,2)
plot(t,v(1:i),'.',t,abs(vs(t)),'linewidth',2);  
title('Speed, m/s');  xlabel('Time, s')
grid; axis tight

Now we add animation to our code:

clear;
g = 9.807;        % gravitational acceleration, m/s^2
A = 1; w = .5;    % amplitude, m, and frequency of surface oscillation, Hz or 1/s 
m = .08; k = .02; % mass of a ball, kg, and coefficient of air resistance, kg/m
H = 2.5;          % initial height of a ball, m
v0 = 0;           % initial speed of a ball, m/s
J = 5;            % maximal number of bounces
T = 5;            % maximal duration of a process, s
dt = 1e-5;        % time integration step, s 
epsh = .0001;     % a threshold height of a ball over surface, m
epsv = .0005;     % a threshold for speed of a ball, m/s
n = round(T/dt);  % maximal number of steps
s = @(t) A*sin(2*pi*w*t);         % surface oscillation function
vs = @(t) A*2*pi*w*cos(2*pi*w*t); % surface velocity function
i = 1;  % steps counter
j = 0;  % bounces counter
h(n) = 0; v(n) = 0;   
h(1) = H; v(1) = v0;  % initial values
bounceindex(J,2)=0; % memory allocation
D=400; % To speed up or slow down the animation, you need
       % to increase or decrease D. With a suitable choice 
       % of D, you can observe the process in real time
while j <= J  
  while i < n % down
    v(i + 1) = (g - k*v(i)^2/m)*dt + v(i);
    h(i + 1) = -v(i)*dt + h(i);
    S = s(dt*(i + 1));    
    if abs(S - h(i + 1)) < epsh % a ball meets a surface
      h(i + 1) = S + epsh;
      v(i + 1) = v(i) + vs(i*dt);
      i = i + 1;
      j = j + 1;
      bounceindex(j,1:2)=[i,j];
      break
    end   
    i = i + 1;
  end 
  
  if j >= J || i >= n
    h(i + 1) = S;
    break 
  end

  while i < n  % up
    v(i + 1) = -(g + k*v(i)^2/m)*dt + v(i);    
    h(i + 1) = v(i)*dt + h(i);
    S = s(dt*(i + 1));  
    if abs(v(i + 1)) < epsv % speed of a ball near zero
      break
    end
     if (abs(S - h(i + 1))) < epsh % a ball meets a surface
       h(i + 1) = S + epsh;
       v(i + 1) = v(i) + vs(i*dt);  
       j = j + 1;
       bounceindex(j,1:2)=[i,j];
       if j >= J || i >= n
         h(i + 1) = S;
         break
       end
     end  
    i = i + 1;
  end
end

figure('menubar', 'none','Name','Bouncing Ball via Euler''s Method')
subplot(2,1,1)
t = dt*(0:i - 1);
plot(t,h(1:i),t,s(dt*(1:i)),'linewidth',2);
title('Height, m');  xlabel('Time, s')
legend('Ball','Surface')
grid; axis tight; 
subplot(2,1,2)
plot(t,v(1:i),'.',t,abs(vs(t)),'linewidth',2);  
title('Speed, m/s');  xlabel('Time, s')
grid; axis tight

Name = 'Animation of Bouncing Ball';
figure('menubar', 'none', 'Name', Name);
ball = rectangle('Position',[4 h(1) 2 2],'Curvature',[1,1],'FaceColor','r');
ground = rectangle('Position',[0 s(t(1))-10 10 10],'FaceColor','k');
elevation = text('Position',[8 h(1)+2]);
surface = text('Position',[8 h(1)+1.5]); 
time = text('Position',[0 h(1)+2]);
height = text('Position',[8 h(1)+2.5]);
bounce = text('Position',[0 h(1)+1.5]);
k = 1;	 
axis off
tic
for i = [2:D:numel(t) - 1,numel(t) - 1]
  set(ball,'Position',[4 h(i) 2 2]);
  set(ground,'Position',[0 s(t(i)) - 10 10 10]);
  set(height,'String',['Height, m: ' num2str(round(h(i - 1),3))]);
  set(time,'String',['Time, s: ' num2str(round(t(i),2))]);
  set(surface,'String',['Surface, m: ' num2str(round(s(t(i)),3))]);
  set(elevation,'String',['Elevation, m: '...
    num2str(round(h(i)-s(t(i)),2))]);
  if i >= bounceindex(k,1)-1
    set(bounce,'String',['Bounce: ' num2str(bounceindex(k,2))]);
    k = k + 1;
  end
  hold off;
  axis([0 10 -2 h(1)+3]);
  drawnow
end
dur=num2str(round(toc,2));
set(gcf,'Name',[Name,'. Duration: ',dur,' s'])