The Nature of Code閱讀心得與Python實作：3.10 The Pendulum

這一節要模擬的是擺(pendulum)這個裝置中，構造最簡單、具有理想化性質的單擺(simple pendulum)。

單擺是由一條下端連接著一個擺錘(bob)，上端懸掛在固定的樞軸(pivot)上的擺臂(arm)組合而成，如圖。當單擺靜止時，擺角，也就是擺臂和鉛垂線之間的夾角，大小會是0；當把把靜止的擺錘往旁邊拉高並放開之後，因為重力及擺臂所提供的張力(tension)的作用，它就會依循著一條弧線來回擺動。

在真實世界中，單擺的運動是在3D空間中進行的，但是，因為我們的模擬世界是個2D的世界，所以我們就只考慮在2D空間中的單擺。

因為單擺是懸掛在樞軸上，所以擺錘在受到重力的作用時，並不會直接掉到地上，而會在重力及擺臂所提供的張力共同作用下，來回擺動。這時候，如果想要模擬擺錘的運動，就不能單純的只考慮重力而已，而必須同時考慮擺角的大小，藉此算出擺錘的角加速度以及角速度來進行模擬。

要算出擺錘的角加速度其實不難，簡單的幾何作圖以及三角函數，還有牛頓第二運動定律，就可以辦到。

參考上圖。作用於擺錘的重力F_g，可以被分解為F_gx及F_gy兩個互相垂直的分量。F_gy跟擺臂提供的張力T，大小相等但方向相反，所以擺錘在運動時，才沒有往外飛出去或向樞軸的方向飛過去；F_gx垂直於擺臂，是促使擺錘擺動的作用力。因為擺臂跟擺錘的運動方向互相垂直，所以不管是T或F_gy，都不會影響擺錘的運動；會影響擺錘運動的，就是F_gx。

為了容易辨識，我們把F_gx改名叫做F_p；又因為它是作用在擺錘上，所以畫圖的時候，可以把起始點移到擺錘的中心點上，這樣更能呈現出它的作用對象是擺錘。

那要怎麼算出F_p呢？假設擺角是θ，由圖可以知道，F_p的大小會是

F_g sinθ

除了大小之外，也必須知道方向。因為F_p的作用方向，是朝著把擺錘拉回靜止時候的位置的緣故，所以當單擺位於鉛垂線右手邊時，F_p<0；而當單擺位於鉛垂線左手邊時，F_p>0。另外，當單擺位於鉛垂線右手邊時，擺角θ的值是正的，所以sinθ>0；而當單擺位於鉛垂線左手邊時，擺角θ的值是負的，所以sinθ<0。因此，根據這些條件以及已知的F_p大小，我們可以得到

F_p = -F_g sinθ

若擺錘的質量是m、重力加速度是g，則

F_g = mg

所以，最後可以得到

F_p = -mg sinθ

有了F_p之後，就可以用牛頓第二定律來找出擺錘的角加速度了。不過，因為擺錘是繞著樞軸旋轉，所以我們用的是牛頓第二定律在這種情況下相對應的版本：

τ = Iα

其中，τ是力矩(torque)，是作用在物體上，使物體繞著轉動軸或支點轉動的物理量；I是轉動慣量(moment of inertia)，是衡量讓擺繞著樞軸轉動的困難度的物理量；α則是角加速度。假設擺長，也就是擺臂長度，是r，則τ和I的計算方式為

τ = rF_p = -rmg sinθ

I = mr²

所以

α = - g sinθ / r

利用這條式子，只要知道擺長、擺角和重力加速度，就可以算出擺錘的角加速度。值得注意的是，在式子中，並未出現擺錘的質量；也就是說，擺錘的角加速度並不會受到擺錘質量大小的影響。

在寫程式時，雖然在真實世界中，地球上的重力加速度是9.8 m/s²，不過，畢竟我們只是要在模擬世界中呈現出擺錘的運動效果，所以就跟在第二章處理萬有引力常數時一樣，在模擬時，我們並不需要去管真正的重力加速度值是多少，只要任意設定一個看起來效果不錯的數字就可以了。

假設變數gravity裡頭放的，是我們設定的重力加速度值、擺角是angle、arm_length是擺長。利用前面導出的擺錘角加速度計算式，以及先前處理角速度和角加速度的方法，模擬的程式可以寫成：

angular_acceleration = -gravity*math.sin(angle)/arm_length
angular_velocity += angular_acceleration
angle += angular_velocity

有了擺角之後，因為擺長是已知的，而且樞軸的位置是我們自己選的，所以就可以算出擺錘的位置。假設樞軸的位置在(x₀, y₀)，而擺角和擺長分別是θ和r，從下圖可以很容易就看出來，擺錘的位置會在

(x₀ + r sinθ, y₀ + r cosθ)

所以，如果pivot=pygame.Vector2(x0, y0)是樞軸位置，計算擺錘位置的程式可以這樣寫

x = arm_length*math.sin(angle)
y = arm_length*math.cos(angle)
bob_position = pivot + pygame.Vector2(x, y)

或者更精簡地寫成

x, y = math.sin(angle), math.cos(angle)
bob_position = pivot + arm_length*pygame.Vector2(x, y)

有了這些之後，就可以設計模擬所需的Pendulum類別了。設計完成的Pendulum類別如下：

class Pendulum:
    def __init__(self, x, y, arm_length, angle_deg):
        # 取得顯示畫面
        self.screen = pygame.display.get_surface()
        
        # 樞軸位置，由傳入的x、y數值決定
        self.pivot = pygame.Vector2(x, y)
        
        # 擺長、擺角
        self.arm_length = arm_length
        self.angle_rad = math.radians(angle_deg)
        
        # 擺錘的初始角加速度及初始角速度
        self.angular_acceleration = 0
        self.angular_velocity = 0
        
        # 擺錘位置
        x, y = math.sin(self.angle_rad), math.cos(self.angle_rad)
        self.bob_position = self.pivot + self.arm_length*pygame.Vector2(x, y)
         
        # 阻尼
        self.damping = 0.995
    
    def update(self):
        # 重力加速度
        gravity = 0.4
        
        self.angular_acceleration = -gravity*math.sin(self.angle_rad)/self.arm_length
        self.angular_velocity += self.angular_acceleration
        self.angle_rad += self.angular_velocity

        x, y = math.sin(self.angle_rad), math.cos(self.angle_rad)
        self.bob_position = self.pivot + self.arm_length*pygame.Vector2(x, y)
        
        self.angular_velocity *= self.damping
        
    def show(self):
        pygame.draw.line(self.screen, (0, 0, 0), self.pivot, self.bob_position)
        pygame.draw.circle(self.screen, (0, 0, 0), self.bob_position, 16)

在設計Pendulum類別時，我們加入了阻尼(damping)，來讓單擺擺動幅度越來越小，這樣子模擬出來的效果，會比較像真實世界中單擺的擺動方式。

先前我們分析所得到的結果，都是針對具理想化性質的單擺所做的分析和推導。這也就是說，我們不僅假設擺錘的質量集中在中心點的位置，同時也沒有考慮樞軸的摩擦力、擺臂的重量、空氣阻力、擺臂是不是會彎曲等等問題；這些在真實世界中，都會對單擺的運動造成影響。不過，畢竟在模擬世界中的單擺，不需要鉅細靡遺地去斤斤計較這麼多，只要它擺動的方式，看起來跟在真實世界中的擺動方式一樣就可以了。在真實世界中的單擺，不管設計再精良，擺錘的擺動幅度，一定會越來越小，最後停止擺動。那要怎麼讓我們模擬出來的單擺，也有這樣子的特性呢？

要讓模擬世界中的單擺擺動幅度越來越小，做法很簡單，只要加個阻尼，讓每一幀的角速度都比原來小一點就可以了。程式可以這樣寫：

angular_velocity *= damping

這裡的damping裡頭放的，是一個小於1的數字。這樣子，當每次計算角速度時，都會比原先的數值小一些，進而使得擺錘的擺動幅度越來越小，到最後停止擺動。

最後，需要提到的是，在設計Pendulum類別時，原書是把計算擺錘位置的程式碼，放在show()方法中，而在這裡，則是放在update()方法中。之所以會如此做，是因為在模擬時，我們不見得會想要畫出擺錘。這時，如果依照原書的設計方式，我們不會執行show()，那當然也就不會知道擺錘的位置在哪裡。所以，把計算擺錘位置的程式寫在update()中，讓處理系統狀態的部分和處理呈現畫面部分的功能各自獨立，這樣子才不會糾纏不清，應該會是比較好的做法；事實上，先前在設計Mover()類別時，就是這麼做的。

Example 3.11: Swinging Pendulum

# python version 3.10.9
import math
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Example 3.11: Swinging Pendulum")

WHITE = (255, 255, 255)

WIDTH, HEIGHT = screen_size = (640, 360)
screen = pygame.display.set_mode(screen_size)

FPS = 60
frame_rate = pygame.time.Clock()

# 樞軸位置
x, y = WIDTH//2, 10

# 擺長
arm_length = 225

# 擺角，單位是「度」
angle_deg = 45

pendulum = Pendulum(x, y, arm_length, angle_deg)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()
            
    screen.fill(WHITE)    
    
    pendulum.update()
    pendulum.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)

Exercise 3.15

# python version 3.10.9
import math
import sys

import pygame  # version 2.3.0

  
pygame.init()

pygame.display.set_caption("Exercise 3.15")

WHITE = (255, 255, 255)

WIDTH, HEIGHT = screen_size = (640, 360)
screen = pygame.display.set_mode(screen_size)

FPS = 60
frame_rate = pygame.time.Clock()

# 單擺數量
n = 5 # number of pendulums

# 樞軸位置
x, y = WIDTH//2, 5

# 擺長
arm_length = 50

# 擺角，單位是「度」
angle = 25  # in degree

pendulums = [Pendulum(x, y, arm_length+7*i, angle+5*i) for i in range(n)]

# 把所有單擺的阻尼都設為1，讓它們永遠都動個不停，這樣看起來會更有趣
for pendulum in pendulums:
    pendulum.damping = 1
    
while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()
            
    screen.fill(WHITE)    
    
    pendulums[0].update()
    for i in range(1,n):
        # 把樞軸放在上一個單擺擺錘的位置上
        pendulums[i].pivot = pendulums[i-1].bob_position
        pendulums[i].update()

    for pendulum in pendulums:
        pendulum.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)

Exercise 3.16

正向力的大小為F_gravity cosθ

Exercise 3.17

假設方塊的質量是m，而重力加速度是g，則

F_gravity = mg

若斜面的摩擦係數是μ，因方塊下滑的力道為F_gravity sinθ，而摩擦力為-μF_gravity cosθ，所以合力為F_gravity(sinθ - μcosθ)。換句話說，當sinθ ≤ μcosθ時，方塊是靜止不動的。

class Box:
    def __init__(self, x, y, mass):
        # 取得顯示畫面的大小
        self.screen = pygame.display.get_surface()
        self.width, self.height = self.screen.get_size()
        
        # 讓傳遞進來的數值來決定物體的質量
        self.mass = mass  
        
        # 物體的質量越大，尺寸就會越大
        self.size = 16*self.mass
        
        # box的傾斜角度
        self.angle_deg = 0
        
        # 垂直於斜面，長度為box邊長一半的向量
        u = pygame.Vector2(self.size/2, 0).rotate(self.angle_deg+45)
                
        # box在斜面上的起始位置
        self.position = pygame.Vector2(x, y) - u
        
        # 物件的初始速度、初始加速度
        self.velocity = pygame.Vector2(0, 0)
        self.acceleration = pygame.Vector2(0, 0)
        
        # 設定box所在surface的格式為per-pixel alpha
        self.surface = pygame.Surface((self.size, self.size), pygame.SRCALPHA)        
            
    def apply_force(self, force):
        self.acceleration += force/self.mass
        
    def update(self):
        self.velocity += self.acceleration
        self.position += self.velocity

        self.acceleration *= 0
        
    def show(self):
        rect = pygame.Rect(0, 0, self.size, self.size)
        pygame.draw.rect(self.surface, (0, 0, 0, 50), rect)
        
        rotated_surface = pygame.transform.rotate(self.surface, self.angle_deg)
        rect_new = rotated_surface.get_rect(center=self.position)
        self.screen.blit(rotated_surface, rect_new)
        
    def check_edges(self):
        if self.position.x - self.size/2 < 0:
            self.position.x = self.size/2
            self.velocity *= 0
    
        if self.position.y + self.size/2 >= self.height:
            self.position.y = self.height - self.size/2
            self.velocity.y = 0

                        
class Incline:
    def __init__(self, x, y, angle_deg):
        # 取得顯示畫面的大小
        self.screen = pygame.display.get_surface()
        self.width, self.height = self.screen.get_size()
        
        # 斜面的傾斜角度，單位是「度」
        self.angle_deg = angle_deg
        
        # 斜面頂部位置
        self.top = pygame.Vector2(x, y)

        # 斜面底部位置    
        dy = self.height - y
        dx = dy/math.tan(math.radians(angle_deg))
        self.bottom = pygame.Vector2(x-dx, self.height)
    
    def on_incline(self, position_x):
        # 假設物體已在斜面上，所以只檢查水平位置
        x1, x2 = self.top.x, self.bottom.x
        x_max, x_min = max(x1, x2), min(x1, x2)
        if x_min <= position_x <= x_max:
            return True
        else:
            return False
        
    def show(self):
        pygame.draw.line(self.screen, (0, 0, 0), self.top, self.bottom)


def friction_force(mu, velocity, N=1):
    if velocity.length() > 0:
        v_hat = velocity.normalize()
    else:
        v_hat = pygame.Vector2(0, 0)  
    
    return -mu*N*v_hat

        
# python version 3.10.9
import math
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Exercise 3.17")

WHITE = (255, 255, 255)

width, height = screen_size = 640, 360
screen = pygame.display.set_mode(screen_size)

FPS = 20
frame_rate = pygame.time.Clock()

# 斜面頂部位置
x, y = 600, 150

# 斜面傾斜角度
angle_deg = 45
angle_rad = math.radians(angle_deg)

incline = Incline(x, y, angle_deg)
slide_direction = (incline.bottom - incline.top).normalize()

# box質量
mass = 3

box = Box(x, y, mass)

# 重力加速度
g = 1

# 斜面及水平面的摩擦係數
mu = 0.8

# 方塊是否正在滑下
sliding = True if math.sin(angle_rad)-mu*math.cos(angle_rad)>0 else False

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()
    
    screen.fill(WHITE)    
    
    if incline.on_incline(box.position.x):
        box.angle_deg = incline.angle_deg

        if sliding:
            slide_force = mass*g*math.sin(angle_rad)*slide_direction
        
            normal_force = mass*g*math.cos(angle_rad)
            friction = friction_force(mu, slide_direction, normal_force)
        
            net_force = slide_force - friction
        else:
            # 方塊為靜止狀態，合力為0
            net_force = pygame.Vector2(0, 0)
    else:
        # 讓方塊滑到底部之後往前翻滾
        if box.position.x + box.size/2 >= incline.bottom.x:
            box.angle_deg += 5
        else:
            box.angle_deg = 90
        
        # 在底部滑動時，如果算出來的摩擦力方向反轉，代表方塊已呈靜止狀態，
        # 而摩擦力已由動摩擦力轉變成靜摩擦力，所以合力應為0
        friction = friction_force(mu, box.velocity, mass*g)
        if friction.x < 0:
            box.velocity = pygame.Vector2(0, 0)
            net_force = pygame.Vector2(0, 0)
        else:
            net_force = friction
            
    box.apply_force(net_force)
    box.update()
    box.check_edges()
    box.show()
    
    incline.show()    
        
    pygame.display.update()
    frame_rate.tick(FPS)