一個地球人在臺北: 一維彈性碰撞次數與圓周率

作者：王一哲
日期：2019/8/14

前言

如下圖所示，假設在水平光滑桌面上有兩個木塊，右側紅色木塊的質量為$m_1$、初速度方向向左，左側綠色木塊的質量為$m_2$、原為靜止，且$m_1 = 100^n m_2$，最左側有固定的牆壁。若木塊之間、木塊與牆壁之間的碰撞皆為一維彈性碰撞，則碰撞次數與$n$的關係如下。

裝置示意圖

n	碰撞次數
1	31
2	314
3	3141
4	31415
5	314159

如果想要了解背後的原理請看這部影片 So why do colliding blocks compute pi? [1]，影片作者是3Blue1Brown。接下來我們拿之前寫過的〈一維彈性碰撞〉程式碼來修改，試著畫出這個現象。

程式及模擬結果

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
n = 4                                             # m1 = 100^n m2
num = 0                                           # 撞擊次數
d2, m2, v2, c2 = 0.2, 1.0, 0.0, color.green       # 左側被撞的木塊寬度, 質量, 初速, 顏色
d1, m1, v1, c1 = 0.4, m2*100**n, -1.0, color.red  # 右側撞人的木塊寬度, 質量, 初速, 顏色
d3, c3 = 0.2, color.blue                          # 左側牆壁的寬度, 地板及牆壁, 顏色
xmax, xmin = 2.0, -2.0                            # x 軸範圍
xrange = xmax - xmin                              # 畫面寬度
t, dt = 0.0, 0.0001                               # 時間, 時間間隔, 單位為s

"""
 2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="1 Dimension Collisions and &pi;", width=800, height=400,
               center=vec(0, 0.2*xmax, 0), background=color.black, range=0.8*xmax, autoscale=False)
# 產生地板
floor = box(size=vec(1.6, 0.04, 0.2)*xrange, pos=vec(0.3*xrange, -0.02*xrange, 0), color=c3)
# 產生牆壁
wall = box(size=vec(d3, 5*d3, 0.2*xrange), pos=vec(xmin + 0.5*d3, 2.5*d3, 0), color=c3)
# 產生左側木塊 b2 並設定初速度
b2 = box(size=vec(d2, d2, d2), pos=vec(0, 0.5*d2, 0), color=c2, m=m2, v=vec(v2, 0, 0))
# 產生右側木塊 b1 並設定初速度
b1 = box(size=vec(d1, d1, d1), pos=vec(xmax - 0.5*d1, 0.5*d1, 0), color=c1, m=m1,  v=vec(v1, 0, 0))
# 產生顯示撞擊次數用的標籤
counter = label(pos=vec(0.5*xmax, 0.8*xmax, 0), text="Collisions: 0", space=50,
                height=24, border=4, box=False, font="monospace")
# 產生顯示質量用的標籤
mass = label(pos=vec(-0.5*xmax, 0.8*xmax, 0),
             text="<i>m</i><sub>1</sub> = 100<sup>%d</sup> <i>m</i><sub>2</sub>" % int(n),
             space=50, height=24, border=4, box=False, font="monospace")
# 繪圖部分
val = abs(sqrt(m1)*v1)*1.1
gd = graph(title="sqrt{m<sub>1</sub>}v<sub>1</sub> - sqrt{m<sub>2</sub>}v<sub>2</sub> plot",
           x=0, y=400, width=400, height=400, xmin=-val, xmax=val, ymin=-val, ymax=val,
           xtitle="sqrt{m<sub>1</sub>}v<sub>1</sub>", ytitle="sqrt{m<sub>2</sub>}v<sub>2</sub>")
fig = gcurve(graph=gd, color=c1)

# 自訂函式，一維彈性碰撞速度公式
def collision(m1, m2, v1, v2):
    v1_prime = (m1-m2)/(m1+m2)*v1 + (2*m2)/(m1+m2)*v2
    v2_prime = (2*m1)/(m1+m2)*v1 + (m2-m1)/(m1+m2)*v2
    return v1_prime, v2_prime

"""
 3. 物體運動部分, 重複執行直到不再發生碰撞為止
"""
while True:
    if(abs(b1.pos.x - b2.pos.x) <= 0.51*(d1 + d2)):
        rate(2000)
        dt = 0.000005
    else:
        rate(1000)
        dt = 0.0005
# 計算木塊間的距離, 若發生碰撞則計算撞後速度
    if(abs(b1.pos.x - b2.pos.x) <= 0.5*(d1 + d2) and b1.v.x < b2.v.x):
        b1.v.x, b2.v.x = collision(b1.m, b2.m, b1.v.x, b2.v.x)
        num += 1
        counter.text = "Collisions: %d" % int(num)
# 計算左側木塊與牆壁的距離, 若發生碰撞則撞後速度反向
    if(abs(b2.pos.x - wall.pos.x) <= 0.5*(d2 + d3) and b2.v.x < 0):
        b2.v.x = -b2.v.x
        num += 1
        counter.text = "Collisions: %d" % int(num)
# 繪圖
    fig.plot(pos=(sqrt(m1)*b1.v.x, sqrt(m2)*b2.v.x))
# 條件成立時停止迴圈
    if(b1.v.x > 0 and b2.v.x > 0 and b1.v.x > b2.v.x): break
# 更新木塊的位置
    b1.pos += b1.v * dt
    b2.pos += b2.v * dt
# 更新時間
    t += dt