Exercise_09 : The Billiard Problem

python physics 陈佩卓 2014301020054

Abstract

Here we consider the problem of a ball moving without friction on a horizontal table. We imagine that there are walls at the edges of the table that reflect the ball perfectly and that there is no friction force between the ball and the table. We can think of this as a billiard ball that moves without friction on a perfect billiard table. The ball is given some initial velocity,and the problem is to calculate and understand the resulting trajectory. This is known as the stadium billiard problem.

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 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.

Main Body

Problem3.31

Study the behavior for other types of tables. One interesting possibility is a square table with a circular interior wall located either in the center, or slightly off-center. Another possibility is an elliptical table.

For the trajectory of the billiard ball in the case of square boundary

Code

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()

When N=2000

For the trajectory of the billiard ball in the case of elliptical boundary

Code

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

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

Code

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

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

Code

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

For the trajectory of the billiard ball in the case of circular boundary.

Code

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

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.

Acknowledgment

Xiongpeiyu