YouTube
Simulation
Model of a Bouncing Ball
Impulse-momentum theorem \begin{eqnarray} \label{eq:im5} m\boldsymbol{v}'=m\boldsymbol{v}+\boldsymbol{I}=m\boldsymbol{v}+I\boldsymbol{n} \end{eqnarray} Collision \begin{eqnarray} \left|\boldsymbol{v}'\cdot\boldsymbol{n}\right|=e\left|\boldsymbol{v}\cdot\boldsymbol{n}\right| \end{eqnarray} Calculating $I$ \begin{eqnarray} I &=& \left| \left(m\boldsymbol{v}'\right) \cdot \boldsymbol{n} \right| + \left| \left(m\boldsymbol{v}\right) \cdot \boldsymbol{n} \right| \nonumber\\ &=& m\left(\left|\boldsymbol{v}' \cdot \boldsymbol{n} \right| + \left| \boldsymbol{v} \cdot \boldsymbol{n} \right|\right) \nonumber\\ &=& m\left(e\left|\boldsymbol{v}\cdot\boldsymbol{n}\right|+\left|\boldsymbol{v}\cdot\boldsymbol{n}\right|\right) \nonumber\\ &=& -m\left(e+1\right) \left( \boldsymbol{v} \cdot \boldsymbol{n} \right) \label{eq:I} \end{eqnarray} Velocity after the collision (inserting (\ref{eq:I}) into (\ref{eq:im5})) \begin{eqnarray} \boldsymbol{v}'=\boldsymbol{v}-(e+1)\left(\boldsymbol{v}\cdot\boldsymbol{n}\right)\boldsymbol{n} \end{eqnarray}
Case 1
\[ f(x, y) = 15^2-\left(x^2+y^2 \right) \]PLT file
# Setting --------------------
reset
set term gif animate delay 5 size 1280, 720
set output "fxy01.gif"
set nokey
set sample 10000
set xr[-20:20]
set yr[-20:20]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics font 'TimesNewRoman, 20'
set xtics 5
set ytics 5
set size ratio -1 # ratio = yrange / xrange
set grid
N = 5
array x[N]
array y[N]
array vx[N]
array vy[N]
array c[N] = ["royalblue", "red", "orange", "green", "black"]
# Parameter --------------------
r = 0.3 # Radius of the ball
g = 9.8 # Gravitational acceleration
dis = 500 # Start to disappear
e = 1.
dt = 0.001 # Time step
dh = dt/6
cut = 80 # Decimation
ep = 1e-4 # for collision detection
time1 = 2 # Stop time4
lim1 = time1/dt/cut
time2 = 40 # Time limit40
lim2 = time2/dt
# Functions --------------------
# Wall
a = 1
b = 1
R = 15
f(x, y) = (x/a)**2 + (y/b)**2
G(x, y) = R**2 - f(x, y)
# Partial derivative of G(x, y)
Gx(x, y) = -(f(x-2*dt, y)-8*f(x-dt, y)+8*f(x+dt, y)-f(x+2*dt, y)) / (12*dt)
Gy(x, y) = -(f(x, y-2*dt)-8*f(x, y-dt)+8*f(x, y+dt)-f(x, y+2*dt)) / (12*dt)
# n = (fx, fy) / sqrt(fx^2+fy^2)
d(x, y) = sqrt(x**2+y**2)
nx(x, y) = Gx(x)/d(Gx(x),Gy(y))
ny(x, y) = Gy(y)/d(Gx(x),Gy(y))
# Inner product of v and n
vn(x,y,vx,vy) = vx*nx(x,y)+vy*ny(x,y)
I(x,y,vx,vy) = -(1+e)*vn(x,y,vx,vy)
# Equations of Motion
rk1(x,y,vx,vy) = vx # dx/dt
rk2(x,y,vx,vy) = vy # dy/dt
rk3(x,y,vx,vy) = 0 # dvx/dt
rk4(x,y,vx,vy) = -g # dvy/dt
# 4th order Runge-Kutta (Define RK_i(x, y, vx, vy))
do for[i=1:4]{
RKi = "RK"
rki = "rk".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, rki, rki, rki, rki)
eval RKi
}
# Elastic
e(e) = sprintf("{/TimesNewRoman:Italic e} = %.2f", e)
# Time
Time(t) = sprintf("{/TimesNewRoman:Italic t} = %.1f s", t)
# Plot --------------------
# Initiate value
t = 0.0 # Time
do for[j=1:N]{
tmpx = 0.8*R*(2*rand(0)-1)
tmpy = 0.8*R*(2*rand(0)-1)
while(f(tmpx, tmpy)>(0.8*R)**2){
tmpx = a*(2*rand(0)-1)
tmpy = b*(2*rand(0)-1)
}
x[j] = tmpx # Position
y[j] = tmpy
vx[j] = 5*(2*rand(0)-1) # Velocity
vy[j] = 5*(2*rand(0)-1)
}
# Draw initial state for lim1 steps
do for[i=1:lim1]{
# Time
set label 1 Time(t) font 'TimesNewRoman, 24' at graph 0.03, 0.04 left
set label 2 e(e) font 'TimesNewRoman, 24' at graph 0.023, 0.09 left
# Wall
set obj 1 circ at 0, 0 size R fs transparent border lc rgb 'black'
# Ball
do for[j=1:N]{
set object j+1 circle at x[j], y[j] size r fc rgb c[j] fs solid front
}
plot 1/0
}
# Update
do for[i=1:lim2]{
# Time
t = t + dt
set label 1 Time(t)
do for[j=1:N]{
# 4th order Runge-Kutta
tmp_x = x[j] + RK1(x[j], y[j], vx[j], vy[j])
tmp_y = y[j] + RK2(x[j], y[j], vx[j], vy[j])
tmp_vx = vx[j] + RK3(x[j], y[j], vx[j], vy[j])
tmp_vy = vy[j] + RK4(x[j], y[j], vx[j], vy[j])
tmp_G = G(tmp_x, tmp_y)
tmp1_G = G(x[j], y[j])
# Collision detection
if(G(x[j], y[j]) >= 0 && G(tmp_x, tmp_y) < 0){
# vec = x_i+1 - x_i
vec_x = tmp_x - x[j]
vec_y = tmp_y - y[j]
nor_v = d(vec_x, vec_y)
while(G(x[j], y[j]) >= 0){
x[j] = x[j] + vec_x / nor_v * ep
y[j] = y[j] + vec_y / nor_v * ep
}
x[j] = x[j] - vec_x / nor_v * ep
y[j] = y[j] - vec_y / nor_v * ep
size = d(tmp_x-x[j], tmp_y-y[j])
tmp_x = x[j] - (1+e)*size* vn(x[j],y[j],vec_x,vec_y) / nor_v * nx(x[j], y[j])
tmp_y = y[j] - (1+e)*size* vn(x[j],y[j],vec_x,vec_y) / nor_v * ny(x[j], y[j])
# Collision
tmp_vx = vx[j] + I(x[j], y[j], vx[j], vy[j]) * nx(x[j], y[j])
tmp_vy = vy[j] + I(x[j], y[j], vx[j], vy[j]) * ny(x[j], y[j])
}
x[j] = tmp_x
y[j] = tmp_y
vx[j] = tmp_vx
vy[j] = tmp_vy
set object 1+N*i+j circ at x[j], y[j] size r fc rgb c[j] fs solid front
set object 1+N*(i-1)+j circ size r*0.01
# Start to disappear
if(i>=dis){
unset object 1+N*(i-dis)+j
}
}
# Decimate and draw
if(i%cut==0){
replot
}
}
set out
GIF file ①
$e=1.00$
Case 2
\[ f(x, y)=15-\left(\left|x\right|+\left|y\right|\right) \]GIF file ②
$e=1.00$
Case 3
\begin{eqnarray*} g(x, y) &=& \left(x^2+y^2-1\right)^3+27x^2y^2\\ f(x, y) &=& g\left(\frac{x}{15}, \frac{y}{15}\right) \end{eqnarray*}GIF file ③
$e=1.00$
Bloopers (in YouTube)
Fall from non-differentiable points
GIF file ④
GIF file ⑤
