Simulation of Bouncing Balls | #3: Bouncing on Explicit Functions (1) [gnuplot]

YouTube

Simulation [gnuplot]

Model of a Bouncing Ball

Impulse-momentum theorem $$\begin{eqnarray} \label{eq:im3} 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:im}$)) $$\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)=y-f(x) \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)=\left(\frac{x}{7}\right)^2-13$

PLT file ①

# Setting --------------------
reset
set term gif animate delay 5 size 1280, 720
set output"ball_fx.gif"
 
set nokey
set sample 10000
set xr[-35:35]
set yr[-15:15]
set xl "{/TimesNewRoman:Italic=24 x}"
set yl "{/TimesNewRoman:Italic=24 y}"
set tics font 'Times New Roman,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     = 0.95            # Elasticity coefficient
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
g(x) = x**2
f(x) = g(x/7.)-13.
W(x, y) = y-f(x)
 
# Partial derivative of f(x)
Wx(x) = -(f(x-2*dt)-8*f(x-dt)+8*f(x+dt)-f(x+2*dt)) / (12*dt) # 4th order central difference formula
Wy(y) = 1
                    
# n = (fx, fy) / sqrt(fx^2+fy^2)
d(x, y)  = sqrt(x**2+y**2)
nx(x, y) = Wx(x)/d(Wx(x), Wy(y))
ny(x, y) = Wy(y)/d(Wx(x), Wy(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
}
 
# Label of displaying parameters
# 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]{
    x[j]  = 13.*(j-3)-2             # Position
    y[j]  = 10.
    vx[j] = 0.                      # Velocity
    vy[j] = 0. 
}
    
# Draw initial state for lim1 steps
do for[i=1:lim1]{
    # Time
    set label 1 Time(t) font 'Times New Roman, 24' at graph 0.03, 0.03
    set label 2 e(e) font 'Times New Roman, 24' at graph 0.023, 0.08
    
    # 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 f(x) lc rgb 'black'
}
 
 
 
# Update
do for[i=1:lim2]{
    # Time
    t = t + dt
    set label 1 Time(t)
    
    # Balls
    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])
        
        # Collision detection
        if(W(x[j], y[j]) >= 0 && W(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)
            
            # Find collision point
            while(W(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
            
            # Calculate bounce direction
            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])
        }
        
        # Update position and velocity of a ball
        x[j]  = tmp_x
        y[j]  = tmp_y
        vx[j] = tmp_vx
        vy[j] = tmp_vy
        
        # Draw a ball
        set object N*i+j circ at x[j], y[j] size r fc rgb c[j] fs solid front
        set object N*(i-1)+j circ size r*0.01    # Old ball turns smaller
 
        # Start to disappear
        if(i>=dis){
            unset object N*(i-dis)+j
        }
    }
    
    # Decimate and draw
    if(i%cut==0){
        replot
    }
}
 
set out 

GIF file ①

$e=0.95$

Case 2: $f(x)=\left|x \right|-5$

GIF file ②

$e=1.00$

Case 3: $f(x) = \frac{x^2}{40}+\cos x-10$

GIF file ③

$e=0.95$

Case 4: $f(x)= \frac{3}{4}x\sin\left(\frac{x^2}{10}+\frac{2}{5}\pi\right)$

GIF file ④

$e=0.95$

Extra

Fall from a non-differentiable point $(f(x)=\left|x \right|-5)$

e=0.50

e=0.80