The Nature of Code閱讀心得與Python實作：2.6 The n-Body Problem

到目前為止，我們所模擬的萬有引力，是一個物體吸引另一個物體，或者是一個物體吸引多個物體。然而，在真實世界中，每個物體都會互相吸引，所以在這一節中，就來把模擬的情境，擴展成多個物體互相吸引。

描述多個物體因萬有引力彼此互相吸引時的運動情形，叫做「多體問題」(n-body problem)。多體問題中，只考慮兩個物體的「雙體問題」(two-body problem)，存有解析解，也就是能夠用數學式子精確地描述這兩個物體的運動情形。然而，只比兩個物體再多一個物體的「三體問題」(three-body problem)，看起來雖然很單純，但它們的運動情形，卻複雜到不存在解析解。劉慈欣著名的科幻小說《三體》，就是以三體問題為引子所展開的故事。

雖然多體問題很複雜，不過利用這一章所提到的方法，倒是可以很粗略地模擬出多體問題的效果。接下來就來設計Body類別。

Body類別的設計很簡單，就是把把Attractor類別中的attract()方法，給複製到Mover類別中，然後稍微修改一下就可以了。修改後的attract()方法，程式如下：

def attract(self, body):
    # 引力方向
    d = self.position - body.position
    r_hat = d.normalize() if d.length() > 0 else pygame.Vector2(0, 0)
            
    distance = d.magnitude()
    # 限制距離數值的範圍，避免極端狀況發生
    if distance < 5:
        distance = 5
    elif distance > 25:
        distance = 25
    
    # 引力大小    
    strength = self.G*self.mass*body.mass/distance**2
    
    # 引力
    force =  strength*r_hat

    # 將引力作用到物體上
    body.apply_force(force) 

這裡要注意一下，這裡的attract()寫法，跟上一節在萬有引力部分最後所採用的attract()寫法不太一樣。上一節的寫法，attract()會傳回計算出來的引力，然後再使用mover物件的apply_force()方法，來將引力作用在mover上；但是這裡的寫法，是直接在attract()內，就將引力作用在要作用的物件上。兩種寫法都能達到目的，不過如果以程式閱讀的角度來看，現在這種寫法，讀起來會比較順暢，接近一般的文字閱讀；通常當我們說「物體A吸引物體B」時，就認定引力已經作用在物體B上了。

下面這個例子，就是利用Body類別建立兩個物件，然後讓它們互相吸引。對於「多體問題」而言，物體最後所呈現出的運動情形，會完全取決於它們的初始條件。所以在這個例子中，分別設定了兩個物體方向相反的初始速度，最後兩個物體會呈現出繞圈圈的運動模式。

Example 2.8: Two-Body Attraction

class Body:
    def __init__(self, x, y, mass, G=1):
        # 取得顯示畫面的大小
        self.screen = pygame.display.get_surface()
        self.width, self.height = self.screen.get_size()
        
        self.G = G
        
        # 讓傳遞進來的數值來決定物體的質量
        self.mass = mass  
        
        # 物體的質量越大，尺寸就會越大
        self.size = (16*self.mass)**0.5
        self.radius = self.size/2
        
        # 讓傳遞進來的數值來決定物體的起始位置
        self.position = pygame.Vector2(x, y)
        
        # 物件的初始速度、初始加速度
        self.velocity = pygame.Vector2(0, 0)
        self.acceleration = pygame.Vector2(0, 0)
        
        # 設定body所在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):
        # 使用具透明度的白色把body所在的surface清空
        self.surface.fill((255, 255, 255, 0))
        
        # 畫出具有透明度的body
        center = pygame.Vector2(self.radius, self.radius)
        pygame.draw.circle(self.surface, (0, 0, 0, 50), center, self.radius)
                
        # 把body所在的surface貼到最後要顯示的畫面上
        self.screen.blit(self.surface, self.position-center)
        
    def attract(self, body):
        # 引力方向
        d = self.position - body.position
        r_hat = d.normalize() if d.length() > 0 else pygame.Vector2(0, 0)
                
        distance = d.magnitude()
        # 限制距離數值的範圍，避免極端狀況發生
        if distance < 5:
            distance = 5
        elif distance > 25:
            distance = 25
        
        # 引力大小    
        strength = self.G*self.mass*body.mass/distance**2
        
        # 引力
        force =  strength*r_hat

        # 將引力作用到物體上
        body.apply_force(force)

  
# python version 3.10.9
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Example 2.8: Two-Body Attraction")

WHITE = (255, 255, 255)

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

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

bodyA = Body(320, 100, 30, 0.5)
bodyB = Body(320, 260, 30, 0.5)

# 設定兩個物體方向相反的初始速度
bodyA.velocity = pygame.Vector2(1, 0)
bodyB.velocity = pygame.Vector2(-1, 0)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill(WHITE)    
    
    # A吸引B；B也吸引A
    bodyA.attract(bodyB)
    bodyB.attract(bodyA)
    
    bodyA.update()
    bodyA.show()
    
    bodyB.update()
    bodyB.show()
        
    pygame.display.update()
    frame_rate.tick(FPS) 

Exercise 2.14

加入一個物體，變成三個物體。

# python version 3.10.9
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Exercise 2.14")

WHITE = (255, 255, 255)

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

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

bodyA = Body(320, 40, 32, 0.5)
bodyB = Body(320, 200, 32, 0.5)
bodyC = Body(320, 320, 32, 0.5)

# 設定A、B方向相反的初始速度，C則一開始靜止不動
bodyA.velocity = pygame.Vector2(2, 0)
bodyB.velocity = pygame.Vector2(-2, 0)
bodyC.velocity = pygame.Vector2(0, 0)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill(WHITE)    
    
    # 物體間互相吸引 
    bodyA.attract(bodyB)
    bodyA.attract(bodyC)
    bodyB.attract(bodyA)
    bodyB.attract(bodyC)
    bodyC.attract(bodyA)
    bodyC.attract(bodyB)
    
    bodyA.update()
    bodyA.show()
    
    bodyB.update()
    bodyB.show()
    
    bodyC.update()
    bodyC.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)

這個練習應該用list搭配迴圈來寫會比較省時省力，而且程式也比較簡潔易讀。接下來就來看看，要怎麼用list搭配迴圈來寫包含多個物體的模擬程式。

假設bodys這個list裡頭有很多個Body物件。當bodys中一個叫bodyA的物件被其他bodys中的物件吸引時，程式這樣寫：

for bodyB in bodys:
    # 自己不會吸引自己
    if bodyB is not bodyA:
        bodyB.attract(bodyA)
        
bodyA.update()
bodyA.show()

每一個在bodys裡頭的物件都這麼做一次，程式這樣寫：

for bodyA in bodys:
    for bodyB in bodys:
        # 自己不會吸引自己
        if bodyB is not bodyA:
            bodyB.attract(bodyA)

for body in bodys:
    body.update()
    body.show()

所以這樣子的程式，就可以模擬多個物體彼此吸引的情況。

這段程式，原書的寫法是把bodyA.update()、bodyA.show()放在bodyA的迴圈裡頭，寫成

for bodyA in bodys:
    for bodyB in bodys:
        # 自己不會吸引自己
        if bodyB is not bodyA:
            bodyB.attract(bodyA)

    bodyA.update()
    bodyA.show()

這種寫法，其實並不是那麼的準確。怎麼說呢？我們來分析一下A、B兩個物體互相吸引時的情況，就可以知道了。要模擬A、B兩個物體互相吸引，我們需要先後處理A吸引B，以及B吸引A兩種情況；當然，我們先前就已經知道，這兩個吸引力，大小相同方向相反。假設A和B間的距離是1，依照原書的寫法，我們會先處理B吸引A，也就是

B.attract(A)

當完成後，就馬上更新A的狀態，也就是

A.update()

當A的狀態更新之後，其位置會因為B的引力作用而有所改變，因而A、B間的距離會縮短，不再是1。所以，當接下來在處理A吸引B時，也就是執行

A.attract(B)

時，因為A、B間的距離已經改變了，並不是原來的1，所以計算出來的A作用在B上的引力大小，與先前計算出來的B作用在A上的引力大小，兩者不會一樣大。這結果顯然是有問題的，不過所幸如果A、B間的距離變化不太大，那在模擬時，應該是很難看出有什麼不對勁的。

Example 2.9: n Bodies

# python version 3.10.9
import random
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Example 2.9: n Bodies")

WHITE = (255, 255, 255)

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

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

bodys = [Body(random.randint(0, WIDTH), 
              random.randint(0, HEIGHT), 
              random.uniform(10, 50),
              0.1)
         for i in range(10)]

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill(WHITE)    
    
    for bodyA in bodys:
        for bodyB in bodys:
            # 自己不會吸引自己
            if bodyB is not bodyA:
                bodyB.attract(bodyA)
    
    for body in bodys:
        body.update()
        body.show()
        
    pygame.display.update()
    frame_rate.tick(FPS) 

Exercise 2.15

要把引力變成斥力，只需要在引力前面加個負號就可以了。

在Body類別中新增repel()方法：

def repel(self, body):
    # 引力方向
    d = self.position - body.position
    r_hat = d.normalize() if d.length() > 0 else pygame.Vector2(0, 0)
            
    distance = d.magnitude()
    # 限制距離數值的範圍，避免極端狀況發生
    if distance < 5:
        distance = 5
    elif distance > 25:
        distance = 25
    
    # 引力大小    
    strength = self.G*self.mass*body.mass/distance**2
    
    # 引力加個負號變成斥力
    force =  -strength*r_hat

    # 將斥力作用到物體上
    body.apply_force(force)      

在主程式中，造一個Body物件當作attractor，並把它的位置設成滑鼠游標的位置。

# python version 3.10.9
import random
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Exercise 2.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()

bodys = [Body(random.randint(0, WIDTH), 
              random.randint(0, HEIGHT),
              random.uniform(10, 30),
              0.05)
         for i in range(10)]

# 在滑鼠游標處造一個attractor
x, y = pygame.mouse.get_pos()
attractor = Body(x, y, 25, 0.5)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill(WHITE)    
    
    # 把attractor的位置設定成滑鼠游標的位置
    attractor.position = pygame.Vector2(pygame.mouse.get_pos())
    
    for bodyA in bodys:
        attractor.attract(bodyA)
        for bodyB in bodys:
            if bodyB is not bodyA:
                bodyB.repel(bodyA)
    
    for body in bodys:
        body.update()
        body.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)

Exercise 2.16

在一個引力很大的物體外圍的圓上，放置一些較小的物體。較小物體的初始速度方向，是所在圓位置的切線方向。調整初始速度的大小，就可以讓這些較小的物體，繞著引力很大的物體進行圓周運動。

半徑為r，圓心位於(0, 0)的圓，在圓上的點(r, 0)，其切線方向的單位向量是(0, 1)或(0, -1)。把位置向量(r, 0)旋轉一個角度，即可得到位於圓上一點的座標；至於該點切線方向的單位向量，則可以將(0, 1)或(0, -1)旋轉相同的角度算出。

# python version 3.10.9
import math
import random
import sys

import pygame  # version 2.3.0


pygame.init()

pygame.display.set_caption("Exercise 2.16")

WHITE = (255, 255, 255)

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

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

attractor = Body(WIDTH/2, HEIGHT/2, 500)

# body數量
n = 30

bodys = []
for i in range(n):
    r = random.uniform(50, 200) # body到attractor的距離
    theta = random.uniform(0, 2*math.pi)
    
    # 建立body之後，再設定其初始位置和速度
    body = Body(0, 0, 25, 0.2)
    body.position = pygame.Vector2(WIDTH/2, HEIGHT/2) + pygame.Vector2(r, 0).rotate_rad(theta)
    body.velocity = 15*pygame.Vector2(0, -1).rotate_rad(theta)
    
    bodys.append(body)
    
while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            sys.exit()

    screen.fill(WHITE)    
    
    for bodyA in bodys:
        attractor.attract(bodyA)
        for bodyB in bodys:
            if bodyB is not bodyA:
                bodyB.attract(bodyA)
    
    for body in bodys:
        body.update()
        body.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)