we consider a problem of a ball moving without friction on a perfect billiard table.Between collisions the velocity is constant so we have
After locatingthe collision point,We must next obtain the unit vector normal to the wall at the point of collision,
.
It is then useful to calculate the components of parallel and perpendicular to the wall.These are just
\vec n)
Hence,the velocity after reflection from the wall is
code
#coding:utf-8
import pylab as pl
import numpy as np
import math
class tabel:
def __init__(self, x0=0.2, y0=0, vx0=-0.2, vy0=-0.5, dtstep=0.001, total_time=50):
self.x = [x0]
self.y = [y0]
self.vx = [vx0]
self.vy = [vy0]
self.time = total_time
self.t = [0]
self.dt = dtstep
self.xbound=[-1]
self.y1bound=[0]
self.y2bound=[0]
def run(self):
for i in range(int(self.time / self.dt)):
self.x.append(self.x[i]+self.vx[i]*self.dt)
self.y.append(self.y[i] + self.vy[i] * self.dt)
self.vx.append(self.vx[i])
self.vy.append(self.vy[i])
if (self.x[i+1]**2+self.y[i+1]**2>1):
xtry=(self.x[i]+self.x[i+1])/2
ytry=(self.y[i]+self.y[i+1])/2
cos=xtry/(xtry**2+ytry**2)**0.5
sin=ytry/(xtry**2+ytry**2)**0.5
verticalx=-(self.vx[i]*cos+self.vy[i]*sin)*cos
verticaly=-(self.vx[i]*cos+self.vy[i]*sin)*sin
#parallelx=self.vx*sin**2-sin*cos*self.vy
# parallely=self.vy*cos**2-self.vx*cos*sin
parallelx=self.vx[i]+verticalx
parallely=self.vy[i]+verticaly
self.vx[i+1]=verticalx+parallelx
self.vy[i+1]=verticaly+parallely
#利用向量變換得到反彈后的速度vx和vy
if (xtry**2+ytry**2>1):
self.x[i+1]=xtry
self.y[i+1]=ytry
continue
else:
self.x[i]=xtry
self.x[i]=ytry
continue
def bound(self):
self.xbound= [-1]
self.y1bound = [0]
self.y2bound = [0]
dx = 0.001
for i in range (2000):
self.xbound.append(self.xbound[i]+dx)
self.y1bound.append((1-self.xbound[i+1]**2)**0.5)
self.y2bound.append(-(1-self.xbound[i+1]**2)**0.5)
def show(self):
pl.plot(self.x,self.y,'.',label='tra')
pl.plot(self.xbound,self.y1bound,'--',label='bound')
pl.plot(self.xbound,self.y2bound,'--',label='bound')
pl.xlabel(u'x')
pl.ylabel(u'y')
a=tabel()
a.run()
a.bound()
a.show()
pl.show()
得到圖像如下圖所示
t=50蜒滩，單位圓下反彈
t=100诗舰，單位圓下反彈
phase space plot
考慮到在兩半圓間增加寬為2α的矩形她肯，并得到動畫，代碼如下
#coding:utf-8
import pylab as pl
import numpy as np
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import animation
#用于制作動畫
class tabel:
def __init__(self, x0=0.2, y0=0, vx0=0.1, vy0=0.3, dtstep=0.01, total_time=300,alpha=0.1):
self.x = [x0]
self.y = [y0]
self.vx = [vx0]
self.vy = [vy0]
self.time = total_time
self.t = [0]
self.dt = dtstep
self.xbound=[-1]
self.y1bound=[0]
self.y2bound=[0]
self.alpha=alpha
def run(self):
for i in range(int(self.time / self.dt)):
self.x.append(self.x[i]+self.vx[i]*self.dt)
self.y.append(self.y[i] + self.vy[i] * self.dt)
self.vx.append(self.vx[i])
self.vy.append(self.vy[i])
self.t.append(self.t[i]+self.dt)
if( self.y[i]>=self.alpha):
if (self.x[i+1]**2+(self.y[i+1]-self.alpha)**2>1):
xtry = self.x[i] + self.vx[i]*self.dt/200
ytry= self.y[i]+self.vy[i]*self.dt/200
dtrebound=self.dt/200
while(xtry**2+(ytry-self.alpha)**2<=1):
xtry = xtry + self.vx[i] * self.dt / 200
ytry = ytry + self.vy[i] * self.dt / 200
dtrebound=dtrebound+self.dt/100
self.t.append(self.t[-1]+dtrebound)
cos=xtry/(xtry**2+(ytry-self.alpha)**2)**0.5
sin=(ytry-self.alpha)/(xtry**2+(ytry-self.alpha)**2)**0.5
verticalx=-(self.vx[i]*cos+self.vy[i]*sin)*cos
verticaly=-(self.vx[i]*cos+self.vy[i]*sin)*sin
parallelx=self.vx[i]+verticalx
parallely=self.vy[i]+verticaly
self.vx[i + 1] = verticalx + parallelx
self.vy[i + 1] = verticaly + parallely#利用向量變換得到反彈后的速度vx和vy
self.x[i+1]=xtry
self.y[i+1]=ytry
elif(self.y[i]<self.alpha and self.y[i]>-self.alpha):
if (abs(self.x[i+1])>1):
xtry = self.x[i] + self.vx[i]*self.dt/100
ytry= self.y[i]+self.vy[i]*self.dt/100
dtrebound=self.dt/100
while(abs(xtry)<=1):
xtry = xtry + self.vx[i] * self.dt / 200
ytry = ytry + self.vy[i] * self.dt / 200
dtrebound=dtrebound+self.dt/100
self.vx[i+1]=-self.vx[i]
self.vy[i+1]=self.vy[i]
elif(self.y[i]<-self.alpha ):
if(self.x[i+1]**2+(self.y[i+1]+self.alpha)**2>1):
xtry = self.x[i] + self.vx[i] * self.dt / 200
ytry = self.y[i] + self.vy[i] * self.dt / 200
dtrebound = self.dt / 100
while (xtry ** 2 + (ytry + self.alpha) ** 2 <=1):
xtry = xtry + self.vx[i] * self.dt / 200
ytry = ytry + self.vy[i] * self.dt / 200
dtrebound = dtrebound + self.dt / 200
cos = xtry / (xtry ** 2 + (ytry+self.alpha) ** 2) ** 0.5
sin = (ytry+self.alpha)/ (xtry ** 2 + (ytry+self.alpha) ** 2) ** 0.5
verticalx = -(self.vx[i] * cos + self.vy[i] * sin) * cos
verticaly = -(self.vx[i] * cos + self.vy[i] * sin) * sin
parallelx = self.vx[i] + verticalx
parallely = self.vy[i] + verticaly
self.vx[i + 1] = verticalx + parallelx
self.vy[i + 1] = verticaly + parallely
self.x[i + 1] = xtry
self.y[i + 1] = ytry
def bound(self):
self.xbound= [-1]
self.y1bound = [self.alpha]
self.y2bound = [-self.alpha]
dx = 0.001
for i in range (2000):
self.xbound.append(self.xbound[i]+dx)
self.y1bound.append(self.alpha+(1-self.xbound[i+1]**2)**0.5)
self.y2bound.append(-self.alpha-(1-self.xbound[i+1]**2)**0.5)
def show(self):
pl.plot(self.x,self.y,'-',label='tra',linewidth=0.1)
pl.plot(self.xbound,self.y1bound,'--')
pl.plot(self.xbound,self.y2bound,'--')
pl.xlabel(u'x')
pl.ylabel(u'y')
pl.ylim(-1.1,1.1)
pl.xlim(-1.1,1.1)
pl.axis('equal')
pl.show()
def drawtrajectory(self):
fig = plt.figure()
ax = plt.axes(title=('Stadium with $\\alpha$ = 0.1, - divergence of two trajectories'),
aspect='equal', autoscale_on=False, xlim=(-1.1, 1.1), ylim=(-1.1, 1.1),
xlabel=('x'),
ylabel=('y'))
line1 = ax.plot([], [], 'b:') # 初始化數(shù)據(jù)次洼，line是軌跡蹭劈，point是軌跡的頭部
point1 = ax.plot([], [], 'bo', markersize=10)
images = []
def init(): # 該函數(shù)用于初始化動畫
line1 = ax.plot([], [], 'b:', markersize=8)
point1 = ax.plot([], [], 'bo', markersize=10)
bound1 = ax.plot([], [], 'k.', markersize=1)
bound2 = ax.plot([], [], 'k.', markersize=1)
return line1, point1, bound1, bound2
def anmi(i): # anmi函數(shù)用于每一幀的數(shù)據(jù)更新酿傍，i是幀數(shù)。
ax.clear()
bound1 = ax.plot(self.xbound, self.y1bound, 'k.', markersize=1)
bound2 = ax.plot(self.xbound, self.y2bound, 'k.', markersize=1)
line1 = ax.plot(self.x[0:20* i], self.y[0:20 * i], 'b:', markersize=8)
point1 = ax.plot(self.x[20 * i - 1:20 * i], self.y[20 * i - 1:20 * i], 'bo', markersize=10)
return line1, point1, bound1, bound2
anmi = animation.FuncAnimation(fig, anmi, init_func=init, frames=500, interval=1, blit=False, repeat=False, )
plt.show()
a=tabel()
a.run()
a.bound()
a.show
a.drawtrajectory()
得到下圖
α=0.1
得到動畫如圖
動畫.gif
