YouTube
Simulation [gnuplot]
Model of a Bouncing Ball
Equations of motion \begin{eqnarray} \label{eq:eom1} \left\{ \begin{array}{l} \frac{dX}{dt}&=&V_X\\ \frac{dY}{dt}&=&V_Y\\ \frac{d{V}_{X}}{dt}&=&-g\sin\alpha\\ \frac{d{V}_{Y}}{dt}&=&-g\cos\alpha \end{array} \right. \end{eqnarray} Inelastic collision \begin{equation} V_{Y}^{\rm{out}}=-\color{red}{e}V_{Y}^{\rm{in}} \end{equation} Rotation matrix \begin{equation} R(\alpha) = \left[ \begin{array}{cc} \cos\alpha & -\sin\alpha \\ \sin\alpha & \cos\alpha \end{array} \right] \end{equation} Rotation using the matrix \begin{eqnarray} \left[\begin{array}{c} x \\ y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}X \\ Y \end{array}\right] &=& \left[ \begin{array}{c} X\cos\alpha - Y\sin\alpha \\ X\sin\alpha + Y\cos\alpha \end{array} \right]\\ \left[\begin{array}{c} v_x \\ v_y\end{array}\right] &=& R(\alpha) \left[\begin{array}{c}V_X \\ V_Y \end{array}\right] &=& \left[ \begin{array}{c} V_X\cos\alpha - V_Y\sin\alpha \\ V_X\sin\alpha + V_Y\cos\alpha \end{array} \right]\\ \end{eqnarray}
BallBouncing.plt
# Setting ----------------------------------------
reset
set term gif animate delay 5 size 1280, 720
set output "BallBouncing.gif"
set nokey
set grid
L = 250 # -L<=x<=L
A = -40 # A<=y<=B
B = 160
set xr[-L:L]
set yr[A:B]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics font 'Times New Roman,20'
set xtics 50
set ytics 40
set size ratio -1
# Parameter ----------------------------------------
g = 9.8 # gravitational acceleration [m/s2]
V = 30 # velocity [m/s]
r = 4.0 # Radius of the ball
a = 5*pi/180 # Slope angle [rad]
b = -85*pi/180 # Projection angle [rad]
dt = 0.001 # Step [s]
dh = dt/6.0
cut = 150 # Decimation
dis = 800 # Start to disappear
N = 5 # The number of balls
array e[N] # The coefficient of restitution
array X[N] # Position of balls
array Y[N]
array VX[N] # Velocity of balls
array VY[N]
array x[N] # Position of balls after a rotation
array y[N]
array vx[N] # Velocity of balls after a rotation
array vy[N]
array color[N] = ["royalblue", "red", "orange", "green", "black"]
cnt = 0 # The number of balls being framed out
# Rotation matrix R(a)
# R[1]=R11, R[2]=R12, R[3]=R21, R[4]=R22
array R[4] = [cos(a), -sin(a), sin(a), cos(a)]
# Functions ----------------------------------------
# Equations of Motion
f1(x, y, vx, vy) = vx # dx/dt
f2(x, y, vx, vy) = vy # dy/dt
f3(x, y, vx, vy) = -g*sin(a) # dvx/dt
f4(x, y, vx, vy) = -g*cos(a) # dvy/dt
# 4th order Runge-Kutta (Define rk_i(x, y, vx, vy))
do for[i=1:4]{
rki = "rk"
fi = "f".sprintf("%d", i)
rki = rki.sprintf("%d(x, y, vx, vy) = (\
k1 = %s(x, y, vx, vy),\
k2 = %s(x + dt*k1/2., y + dt*k1/2., vx + dt*k1/2., vy + dt*k1/2.),\
k3 = %s(x + dt*k2/2., y + dt*k2/2., vx + dt*k2/2., vy + dt*k2/2.),\
k4 = %s(x + dt*k3, y + dt*k3, vx + dt*k3, vy + dt*k3),\
dh * (k1 + 2*k2 + 2*k3 + k4))", i, fi, fi, fi, fi)
eval rki
}
# (X,Y) ->(x,y) (x=R(a)X)
x(X, Y) = R[1]*X + R[2]*Y
y(X, Y) = R[3]*X + R[4]*Y
# (x,y) -> (X,Y) (X=R(-a)x)
X(x, y) = R[4]x - R[2]*y
Y(x, y) = -R[3]*x + R[1]*y
# Time
Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f s", t)
# Plot ----------------------------------------
# Initial Value
t = 0
b_inc = 3
do for[i=1:N]{
X[i] = X(150, 60)
Y[i] = Y(150, 60)
VX[i] = V*cos(b + (b_inc*(i-int(N/2))) * pi/180)
VY[i] = V*sin(b + (b_inc*(i-int(N/2))) * pi/180)
e[i] = 0.3+(1.0-0.3)/N*i # float division
# Vector rotation (Rotation matrix R(a))
x[j] = x(X[j] , Y[j] )
y[j] = y(X[j] , Y[j] )
vx[j] = x(VX[j], VY[j])
vy[j] = y(VX[j], VY[j])
}
# Draw initiate state for 70 steps
do for [i=1:70] {
# Time
set label 1 Time(t) left at graph 0.01, 0.93 font 'TimesNewRoman, 25'
# Balls
do for[j=1:N]{
set obj j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
}
# Draw ground and balls
plot tan(a)*x with filledcurve x1 lc rgb 'gray50'
}
# Update until all of balls are framed out
i = 0 # The number of loops
while(1){
i = i +1
t = t + dt
set label 1 Time(t)
do for[j=1:N]{
# 4th order Runge-Kutta
temp_X = X[j] + rk1(X[j], Y[j], VX[j], VY[j])
temp_Y = Y[j] + rk2(X[j], Y[j], VX[j], VY[j])
temp_VX = VX[j] + rk3(X[j], Y[j], VX[j], VY[j])
temp_VY = VY[j] + rk4(X[j], Y[j], VX[j], VY[j])
X[j]=temp_X; Y[j]=temp_Y; VX[j]=temp_VX; VY[j]=temp_VY
# Judge whether balls bounce or not
if(Y[j]<r){
Y[j] = r
VY[j] = -e[j]*VY[j] # inelastic collision
}
# Vector rotation (Rotation matrix R(a))
x[j] = x(X[j] , Y[j] )
y[j] = y(X[j] , Y[j] )
vx[j] = x(VX[j], VY[j])
vy[j] = y(VX[j], VY[j])
# Update objects
set obj N*i+j circ at x[j], y[j] size r fc rgb color[j] fs solid noborder
# Make old objects trajectory of the ball
set obj N*(i-1)+j at x[j], y[j] size 0.1
# Start to disappear
if(i>=dis){
unset object N*(i-dis)+j
}
}
# Decimate and plot
if(i%cut==0){
replot
}
# Count the number of balls being framed out
do for[j=1:N]{
if(x[j]<-L-r){
cnt = cnt + 1
}
}
# Exit from the loop when all balls being framed out
if(cnt == N){
break
} else {
cnt =0
}
}
set out
BallBouncing.gif
$V=30\mathrm{m/s},\ a=5^\circ,\ b=-85^\circ,\ e_i=0.3+0.14i,\ b_{\mathrm{inc}}=3$$t=0: X_j=200\cos\left(a\right)+60\sin\left(a\right),\ Y_{j}=-200\sin\left(a\right)+60\cos\left(a\right)\\ \left(t=0: x_{j}=200,\ y_{j}=60\right)$
BallBounding2.gif
$V=30\mathrm{m/s},\ a=15^\circ,\ b=-75^\circ,\ e=0.7, b_{\mathrm{inc}}=3$$t=0: X_j=200\cos\left(a\right)+120\sin\left(a\right),\ Y_{j}=-200\sin\left(a\right)+120\cos\left(a\right)\\ \left(t=0: x_{j}=200,\ y_{j}=120\right)$
BallBounding3.gif
$V=30\mathrm{m/s},\ a=5^\circ,\ b=15^\circ,\ e=0.8, b_{\mathrm{inc}}=9$$t=0: X_j=-0.8L,\ Y_{j}=0.25B$
Extra
BallBounding4.gif
$V=30\mathrm{m/s},\ a=5^\circ,\ b=-75^\circ,\ e=0.8, b_{\mathrm{inc}}=5$$t=0: X_j=60\sin\left(a\right),\ Y_{j}=60\cos\left(a\right)\\ \left(t=0: x_{j}=0,\ y_{j}=60\right)$
