Exercise_09 : The Billiard Problem

Problem 3.31

冉峰 弘毅班 2014301020064

Abstract

• The billiard system can also be a chaotic system. I will try to solve the trajectories of different kinds of tables and show the phase diagram to see whether the system is chaotic or not. As the previous dicusssion of two kinds of situation in the chaos, I will consider one more chaotice model in this paper. Here I consider the problem of a ball moving without friction on a horizontal table and dicsuss the phase plot,the trajectory plot ant the divergence plot with the stadium ball.

Background

• In this problem, I will not consider the effect of friction, which means the billiard ball will move without friction on a perfect billiard table.

The moving

• Except for the collisions with the walls, the motion of the billiard is quite simple. Between collisions the velocity is constant so we have ：
  $$\frac{dx}{dt} = v_x$$
  $$\frac{dy}{dt} = v_y$$
  where $v_x$ and $v_y$ change only through collisions with the walls. These equations can be solved using our Euler algorithm.

The collision

• It is then useful to calculate the components of $v_{i}$ parallel and perpendicular to the wall.
  These are just:
  $$\vec{v}_{i,\bot} = (\vec{v}_i\cdot\hat{n})\hat{n}$$
  $$\vec{v}_{i,\|} = \vec{v}_{i} - \vec{v}_{i,\bot}$$
  After reflection, the velocity becomes:
  $$\vec{v}_{f,\bot} = -\vec{v}_{i,\bot}$$
  $$\vec{v}_{f,\|} = \vec{v}_{i,\|}$$

The Chaos Theory

• We know Chaos theory is the field of study in mathematics that studies the behavior of dynamical systems that are highly sensitive to initial conditions—a response popularly referred to as the butterfly effect. Small differences in initial conditions (such as those due to rounding errors in numerical computation) yield widely diverging outcomes for such dynamical systems, rendering long-term prediction impossible in general. This happens even though these systems are deterministic, meaning that their future behavior is fully determined by their initial conditions, with no random elements involved. In other words, the deterministic nature of these systems does not make them predictable. This behavior is known as deterministic chaos, or simply chaos. The theory was summarized by Edward Lorenz as:

The Lorenz Attractor

• Having been aware of them,we can safely say that it will not be hard for us to draw all the plot we needed.One points need to be mentioed.A hall mark of a chaotic system is an extreme sensitivity to intial conditions.This property is also found in the Lorentz attractor:
  $$\frac{dx}{dt}=\sigma(y-x)$$
  $$\frac{dy}{dt}=-xz+rx-y$$
  $$\frac{dz}{dt}=xy-bz$$
  The Lorenz variables $x,y,z$ are derived from the temperature, density and velocity variables in the original Navier-Stokes equations, and the parameters $\sigma,r,b$ are measures of the temperature difference across the fluid and other fluid parameters.

The Main Body

Code for the trajectory of the billiard ball in the case of square boundary

import matplotlib.pyplot as plt
import numpy as np
class billiard_rectangular:
    def __init__(self,x_0,y_0,vx_0,vy_0,N,dt):
        self.x_0 = x_0
        self.y_0 = y_0
        self.vx_0 = vx_0
        self.vy_0 = vy_0
        self.N = N
        self.dt = dt
    def motion_calculate(self):
        self.x = []
        self.y = []
        self.vx = []
        self.vy = []
        self.t = [0]
        self.x.append(self.x_0)
        self.y.append(self.y_0)
        self.vx.append(self.vx_0)
        self.vy.append(self.vy_0)
        for i in range(1,self.N):
            self.x.append(self.x[i - 1] + self.vx[i - 1]*self.dt)
            self.y.append(self.y[i - 1] + self.vy[i - 1]*self.dt)
            self.vx.append(self.vx[i - 1])
            self.vy.append(self.vy[i - 1])
            if (self.x[i] < -1.0):
                self.x[i],self.y[i] = self.correct('x>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.x[i] > 1.0):
                self.x[i],self.y[i] = self.correct('x<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.y[i] < -1.0):
                self.x[i],self.y[i] = self.correct('y>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i]
            elif(self.y[i] > 1.0):
                self.x[i],self.y[i] = self.correct('y<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i] 
            else:
                pass
            self.t.append(self.t[i - 1] + self.dt)
        return self.x, self.y
    def correct(self,condition,x,y,vx,vy):
        vx_c = vx/100.0
        vy_c = vy/100.0
        while eval(condition):
            x = x + vx_c*self.dt
            y = y + vy_c*self.dt
        return x-vx_c*self.dt,y-vy_c*self.dt
    def reflect(self):
        pass
    def plot(self):
        plt.figure(figsize = (8,8))
        plt.xlim(-1,1)
        plt.ylim(-1,1)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.plot(self.x,self.y,'g')
        plt.show()

Result

• We can chose different N for $N=5000,2000,500,200$
  
  
  
  

• Actually if the ratio of $v_x$ and $v_y$ is rational, the square cannot be fully filled, for the quantity of irrational number is larger

Vpython code

from visual import *
side = 4.0
thk = 0.3
s2 = 2*side - thk
s3 = 2*side + thk
wallR = box (pos=( side, 0, 0), size=(thk, s2, s3),  color = color.green)
wallL = box (pos=(-side, 0, 0), size=(thk, s2, s3),  color = color.green)
wallB = box (pos=(0, -side, 0), size=(s3, thk, s3),  color = color.green)
wallT = box (pos=(0,  side, 0), size=(s3, thk, s3),  color = color.green)
wallG = box (pos=(0,0,-2*side),size=(s3,s3,s3),color=color.green)
ball = sphere(pos=(0.2,0,0), color=color.white, radius = 0.5, make_trail=True, retain=200)
ball.trail_object.radius = 0.07
ball.v = vector(0.5,0.3,0.4)
side = side - thk*0.5 - ball.radius
dt = 0.5
t=0.0
while True:
  rate(100)
  t = t + dt
  ball.pos = ball.pos + ball.v*dt
  if not (side > ball.x > -side):
    ball.v.x = -ball.v.x
  if not (side > ball.y > -side):
    ball.v.y = -ball.v.y
  if not (side>ball.z>-side):
      ball.v.z=-ball.v.z

Result

• In 2D
  
• In 3D
  

Code for the trajectory of the billiard ball in the case of circular boundary

class billiard_circle(billiard_rectangular):
    def motion_calculate(self):
        self.x = []
        self.y = []
        self.vx = []
        self.vy = []
        self.t = [0]
        self.x.append(self.x_0)
        self.y.append(self.y_0)
        self.vx.append(self.vx_0)
        self.vy.append(self.vy_0)
        for i in range(1,self.N):
            self.x.append(self.x[i - 1] + self.vx[i - 1]*self.dt)
            self.y.append(self.y[i - 1] + self.vy[i - 1]*self.dt)
            self.vx.append(self.vx[i - 1])
            self.vy.append(self.vy[i - 1])
            if (np.sqrt( self.x[i]**2+self.y[i]**2 ) > 1.0):
                self.x[i],self.y[i] = self.correct('np.sqrt(x**2+y**2) < 1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i],self.vy[i] = self.reflect(self.x[i],self.y[i],self.vx[i - 1], self.vy[i - 1])
            self.t.append(self.t[i - 1] + self.dt)
        return self.x, self.y
    def reflect(self,x,y,vx,vy):
        module = np.sqrt(x**2+y**2)  ### normalization
        x = x/module
        y = y/module
        v = np.sqrt(vx**2+vy**2)
        cos1 = (vx*x+vy*y)/v
        cos2 = (vx*y-vy*x)/v
        vt = -v*cos1
        vc = v*cos2 
        vx_n = vt*x+vc*y
        vy_n = vt*y-vc*x
        return vx_n,vy_n
    def plot(self):
        plt.figure(figsize = (8,8))
        plt.xlim(-1,1)
        plt.ylim(-1,1)
        plt.xlabel('x')
        plt.ylabel('y')
        self.plot_boundary()
        plt.plot(self.x,self.y,'b')
        plt.show()
    def phase_space_plot(self):
        plt.figure(figsize = (16,8))
        plt.subplot(121)
        plt.xlabel('x')
        plt.ylabel(r'$v_x$')
        plt.scatter(self.x,self.vx, s=1)
        plt.subplot(122)
        plt.xlabel('y')
        plt.ylabel(r'$v_x$')
        plt.scatter(self.y,self.vx, s=1)
        plt.show()    
    def phase_plot(self):
        plt.figure(figsize = (8,8))
        record_x = []
        record_vx = []
        for i in range(len(self.x)):
            if (abs(self.y[i] - 0)<0.001):
                record_vx.append(self.vx[i])
                record_x.append(self.x[i])
        plt.xlabel('x')
        plt.ylabel(r'$v_x$')
        plt.scatter(record_x,record_vx,s=1)
        plt.show()
    def plot_boundary(self):
        theta = 0
        x = []
        y = []
        while theta < 2*np.pi:
            x.append(np.cos(theta))
            y.append(np.sin(theta))
            theta+= 0.01
        plt.plot(x,y)

Result

• We can chose different N for $N=5000,2000,500,200$

• It shows that the billiard ball cannot reach the center part, which is determined by the position of tha ball and direction of initial. velocity.

This is the phase diagram of $(x,v_x)$ and $(y,v_y)$

Code for the trajectory of the billiard ball in the case of elliptical boundary

class billiard_ellipse(billiard_circle):   ### x^2/3+y^2/2 = 1
    def motion_calculate(self):
        self.x = []
        self.y = []
        self.vx = []
        self.vy = []
        self.t = [0]
        self.x.append(self.x_0)
        self.y.append(self.y_0)
        self.vx.append(self.vx_0)
        self.vy.append(self.vy_0)
        for i in range(1,self.N):
            self.x.append(self.x[i - 1] + self.vx[i - 1]*self.dt)
            self.y.append(self.y[i - 1] + self.vy[i - 1]*self.dt)
            self.vx.append(self.vx[i - 1])
            self.vy.append(self.vy[i - 1])
            if (self.x[i]**2/3+self.y[i]**2/2 > 1.0):
                self.x[i],self.y[i] = self.correct('x**2/3+y**2/2 < 1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i],self.vy[i] = self.reflect((2./3)*self.x[i],self.y[i],self.vx[i - 1], self.vy[i - 1])
            self.t.append(self.t[i - 1] + self.dt)
        return self.x, self.y       
    def plot(self):
        plt.figure(figsize = (8,6))
        plt.xlim(-2.2,2.2)
        plt.ylim(-2,2)
        plt.xlabel('x')
        plt.ylabel('y')
        self.plot_boundary()
        plt.plot(self.x,self.y,'b')
        plt.show()
    def plot_boundary(self):
        theta = 0
        x = []
        y = []
        while theta < 2*np.pi:
            x.append(np.sqrt(3)*np.cos(theta))
            y.append(np.sqrt(2)*np.sin(theta))
            theta+= 0.01
        plt.title(r'Elliptical stadium  $\frac{x^2}{4}+\frac{y^2}{3} = 3$')
        plt.plot(x,y)

Result

• We can chose different N for $N=5000,2000,500,200$ and chose different ellipticals.
  

• It shows that the billiard ball cannot reach the left and the right part, which is determined by the position of tha ball and direction of initial.

• This is Poincare section by scattering the points only when the billiard ball crosses the y = 0 axis, which is a ellipse-liked curve
  

Code for the trajectory of the billiard ball in the case of square with an inner circle boundary, and the circle is in the center

class billiard_innercircle(billiard_circle):
    def motion_calculate(self):
        self.x = []
        self.y = []
        self.vx = []
        self.vy = []
        self.t = [0]
        self.x.append(self.x_0)
        self.y.append(self.y_0)
        self.vx.append(self.vx_0)
        self.vy.append(self.vy_0)
        for i in range(1,self.N):
            self.x.append(self.x[i - 1] + self.vx[i - 1]*self.dt)
            self.y.append(self.y[i - 1] + self.vy[i - 1]*self.dt)
            self.vx.append(self.vx[i - 1])
            self.vy.append(self.vy[i - 1])
            if (self.x[i] < -1.0):
                self.x[i],self.y[i] = self.correct('x>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.x[i] > 1.0):
                self.x[i],self.y[i] = self.correct('x<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.y[i] < -1.0):
                self.x[i],self.y[i] = self.correct('y>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i]
            elif(self.y[i] > 1.0):
                self.x[i],self.y[i] = self.correct('y<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i] 
            elif(self.x[i]**2+self.y[i]**2<0.01):
                self.x[i],self.y[i] = self.correct('x**2+y**2 > 0.01',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i],self.vy[i] = self.reflect(self.x[i],self.y[i],self.vx[i - 1], self.vy[i - 1])
            self.t.append(self.t[i - 1] + self.dt)
        return self.x, self.y    
    def plot_boundary(self):
        theta = 0
        x = []
        y = []
        while theta < 2*np.pi:
            x.append(0.1*np.cos(theta))
            y.append(0.1*np.sin(theta))
            theta+= 0.01
        plt.plot(x,y)

Result

• We can chose different N for $N=5000,500$
  

• This is Poincare section by scattering the points only when the billiard ball crosses the y = 0 axis.
  

Code for the trajectory of the billiard ball in the case of square with an inner circle boundary and the circle is slightly off-center

class billiard_innercircle2(billiard_circle):
    def motion_calculate(self):
        self.x = []
        self.y = []
        self.vx = []
        self.vy = []
        self.t = [0]
        self.x.append(self.x_0)
        self.y.append(self.y_0)
        self.vx.append(self.vx_0)
        self.vy.append(self.vy_0)
        for i in range(1,self.N):
            self.x.append(self.x[i - 1] + self.vx[i - 1]*self.dt)
            self.y.append(self.y[i - 1] + self.vy[i - 1]*self.dt)
            self.vx.append(self.vx[i - 1])
            self.vy.append(self.vy[i - 1])
            if (self.x[i] < -1.0):
                self.x[i],self.y[i] = self.correct('x>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.x[i] > 1.0):
                self.x[i],self.y[i] = self.correct('x<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i] = - self.vx[i]
            elif(self.y[i] < -1.0):
                self.x[i],self.y[i] = self.correct('y>-1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i]
            elif(self.y[i] > 1.0):
                self.x[i],self.y[i] = self.correct('y<1.0',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vy[i] = -self.vy[i] 
            elif((self.x[i] - 0.05)**2+self.y[i]**2<0.01):
                self.x[i],self.y[i] = self.correct('(x - 0.05)**2+y**2 > 0.01',self.x[i - 1], self.y[i - 1], self.vx[i - 1], self.vy[i - 1])
                self.vx[i],self.vy[i] = self.reflect(self.x[i]-0.05,self.y[i],self.vx[i - 1], self.vy[i - 1])
            self.t.append(self.t[i - 1] + self.dt)
        return self.x, self.y    
    def plot_boundary(self):
        theta = 0
        x = []
        y = []
        plt.title(r'$(x-0.5)^2+y^2 = 0.01$')
        while theta < 2*np.pi:
            x.append(0.1*np.cos(theta)+0.5)
            y.append(0.1*np.sin(theta))
            theta+= 0.01
        plt.plot(x,y)    

Result

• We can chose different N for $N=5000,500$
  

• Which is nearly the same as the former case.

Conclusion

• Ths chaotic system will present when the boundary is composed by some shapes (different or same). If the boundary is a simple regular shape, the chaotic situation cannot be observed.

Thanks For

• The book of Computational Physics, chapter 3
• The Lesson plan Chapter 3 of Cai Hao
• Vpython code from Chen feng