计算物理第九次作业Problem3.30

摘要

上次作业研究了参数对混沌现象的影响，而这次作业着重于研究桌球的碰撞问题。

背景

已知运动方程为

正文

1、方形边界
程序为：

import matplotlib.pyplot as plt
vx = 3
vy = 3.5
x = [0.2]
y = [0]
dt = 0.01
for i in range(10000):
    x.append(x[-1]+vx*dt)
    y.append(y[-1]+vy*dt)
    t = i*dt
    if x[-1]>1:
        r=(x[-2]-1)/(1-x[-1])
        y[-1]=(y[-2]+r*y[-1])/(r+1)
        x[-1]=1
        vx=-vx
    if x[-1]<-1:
        r=-(1+x[-2])/(x[-1]+1)
        y[-1]=(y[-2]+r*y[-1])/(r+1)
        x[-1]=-1
        vx=-vx
    if y[-1]>1:
        r=(y[-2]-1)/(1-y[-1])
        x[-1]=(x[-2]+r*x[-1])/(r+1)
        y[-1]=1
        vy=-vy
    if y[-1]<-1:
        r=-(1+y[-2])/(y[-1]+1)
        x[-1]=(x[-2]+r*x[-1])/(r+1)
        y[-1]=-1
        vy=-vy
plt.plot(x,y)
plt.xlabel('x')
plt.ylabel('y')

运行结果为

改变速度方向，可得到：


2、圆形边界
(1)一般圆桌碰撞
程序为：

import matplotlib.pyplot as plt
import numpy as np
v = 2
a = 0.01
x = [0.2]
y = [0]
vx = 2
vy = 1.6
dt = 0.01
for i in range(10000):
    x.append(x[-1]+vx*dt)
    y.append(y[-1]+vy*dt)  
    t = i*dt
    if (x[-1]*x[-1]+y[-1]*y[-1])>1:
        x[-1] = x[-2] + vx/100*dt
        y[-1] = y[-2] + vy/100*dt
        m = np.sqrt(x[-1]**2+y[-1]**2)  
        x[-1] = x[-1]/m
        y[-1] = y[-1]/m
        v = np.sqrt(vx**2+vy**2)
        cos1 = (vx*x[-1]+vy*y[-1])/v
        cos2 = (vx*y[-1]-vy*x[-1])/v
        vt = -v*cos1
        vc = v*cos2 
        vx = vt*x[-1]+vc*y[-1]
        vy = vt*y[-1]-vc*x[-1]
plt.figure(figsize = (5,5))
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.plot(x,y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Circular stadium - trajectory')

运行结果为：

(2)从x轴切开距离为2ar
程序为：

import matplotlib.pyplot as plt
import numpy as np
v = 2
a = 0.01
x = [0.2]
y = [0]
vx = 2
vy = 1.6
dt = 0.01
for i in range(5000):
    x.append(x[-1]+vx*dt)
    y.append(y[-1]+vy*dt)  
    t = i*dt
    if (x[-1]*x[-1]+y[-1]*y[-1])>1 and y[-1]>0.01:
        while (x[-1]*x[-1]+y[-1]*y[-1])<1:
            x[-1] = x[-2] + vx/100*dt
            y[-1] = y[-2] + vy/100*dt
        m = np.sqrt(x[-1]**2+y[-1]**2)  
        x[-1] = x[-1]/m
        y[-1] = (y[-1]-0.01)/m+0.01
        v = np.sqrt(vx**2+vy**2)
        cos1 = (vx*x[-1]+vy*(y[-1]-0.01))/v
        cos2 = (vx*(y[-1]-0.01)-vy*x[-1])/v
        vt = -v*cos1
        vc = v*cos2 
        vx = vt*x[-1]+vc*(y[-1]-0.01)
        vy = vt*(y[-1]-0.01)-vc*x[-1]
    elif (x[-1]*x[-1]+y[-1]*y[-1])>1 and y[-1]<-0.01:
        while (x[-1]*x[-1]+y[-1]*y[-1])<1:
            x[-1] = x[-2] + vx/100*dt
            y[-1] = y[-2] + vy/100*dt
        m = np.sqrt(x[-1]**2+y[-1]**2)  
        x[-1] = x[-1]/m
        y[-1] = (y[-1]+0.01)/m-0.01
        v = np.sqrt(vx**2+vy**2)
        cos1 = (vx*x[-1]+vy*(y[-1]+0.01))/v
        cos2 = (vx*(y[-1]+0.01)-vy*x[-1])/v
        vt = -v*cos1
        vc = v*cos2 
        vx = vt*x[-1]+vc*(y[-1]+0.01)
        vy = vt*(y[-1]+0.01)-vc*x[-1]
    elif x[-1]<-1 and y[-1]>-0.01 and y[-1]<0.01:
        x[-1] = x[-2] + vx/100*dt
        y[-1] = y[-2] + vy/100*dt
        vx=-vx
    elif x[-1]>1 and y[-1]>-0.01 and y[-1]<0.01:
        x[-1] = x[-2] + vx/100*dt
        y[-1] = y[-2] + vy/100*dt
        vx=-vx
plt.figure(figsize = (5,5))
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.plot(x,y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Stadium billiard  a=0.01')

运行结果为：

结论

对于一般圆桌碰撞，无混沌现象；而当a=0.01时，图像会出现不规则的现象。

致谢

感谢周璟怡学姐的耐心帮助。