Simulation of Bouncing Balls | #5: Bouncing on Implicit Functions (1) [gnuplot]

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}$$
Wall function $$\begin{equation} W(x, y)=-f(x, y) \end{equation}$$ Inequality constraint $$\begin{equation} W(x, y)\geq 0 \end{equation}$$ Unit normal vector $$\begin{equation} \boldsymbol{n}=\frac{\nabla W}{\left|\nabla W\right|} \end{equation}$$ Collision detection $$\begin{equation} \boldsymbol{v}'=\boldsymbol{v}-\left(e+1\right)\left(\boldsymbol{v}\cdot \frac{\nabla W}{\left|\nabla W\right|}\right)\frac{\nabla W}{\left|\nabla W\right|} \end{equation}$$

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 ⑤