The Nature of Code閱讀心得筆記——使用Python實作: 2.6 The n-Body Problem

2.6 The n-Body Problem

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

先前我們所寫的Mover和Attractor兩個類別，在模擬時，就只有Attractor類別的物件會吸引Mover類別的物件。但既然在真實世界中，每個物體都會互相吸引，所以我們不再使用Mover和Attractor這兩個類別，而是設計一個新的Body類別，讓Body類別的不同物件互相吸引。

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

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

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

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)

下面這個例子，就是利用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.radius = (16*self.mass)**0.5
        
        # 讓傳遞進來的數值來決定物體的起始位置
        self.position = pygame.Vector2(x, y)
        
        self.velocity = pygame.Vector2(0, 0)
        self.acceleration = pygame.Vector2(0, 0)
        
        # 設定body所在surface的格式為per-pixel alpha
        self.top_surface = pygame.Surface((self.width, self.height), 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.top_surface.fill((255, 255, 255, 0))
        
        # 畫出具有透明度的body
        pygame.draw.circle(self.top_surface, (0, 0, 0, 50), self.position, self.radius)
                
        # 把body所在的surface貼到最後要顯示的畫面上
        self.screen.blit(self.top_surface, (0, 0))
        
    def check_edges(self):
        if self.position.x > self.width:
            self.position.x = self.width
            self.velocity.x = -self.velocity.x
        elif self.position.x < 0:
            self.position.x = 0
            self.velocity.x = -self.velocity.x
    
        if self.position.y > self.height:
            self.position.y = self.height
            self.velocity.y = -self.velocity.y
   
    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)        
        
    
import sys
import pygame

pygame.init()

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

FPS = 60
WHITE = (255, 255, 255)

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

frame_rate = pygame.time.Clock()

bodyA = Body(320, 40, 8)
bodyB = Body(320, 200, 8)

# 設定兩個物體方向相反的初始速度
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

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

import sys
import pygame

pygame.init()

pygame.display.set_caption("Exercise 2.14")

FPS = 60
WHITE = (255, 255, 255)

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

frame_rate = pygame.time.Clock()

bodyA = Body(320, 40, 8)
bodyB = Body(320, 200, 8)
bodyC = Body(320, 320, 8)

# 設定A、B方向相反的初始速度，C則一開始靜止不動
bodyA.velocity = pygame.Vector2(1, 0)
bodyB.velocity = pygame.Vector2(-1, 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)

    bodyA.update()
    bodyA.show()

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

Example 2.9: n Bodies

import random
import sys
import pygame

pygame.init()

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

FPS = 60
WHITE = (255, 255, 255)

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

frame_rate = pygame.time.Clock()

bodys = [Body(random.randint(0, WIDTH), 
              random.randint(0, HEIGHT), 
              random.uniform(0.1, 5),
              0.4)
         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)
    
        bodyA.update()
        bodyA.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，並把它的位置設成滑鼠游標的位置。

import random
import sys
import pygame

pygame.init()

pygame.display.set_caption("Exercise 2.15")

FPS = 60
WHITE = (255, 255, 255)

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

frame_rate = pygame.time.Clock()

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

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

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)
    
        bodyA.update()
        bodyA.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)旋轉相同的角度算出。

import math
import random
import sys
import pygame

pygame.init()

pygame.display.set_caption("Exercise 2.16")

FPS = 60
WHITE = (255, 255, 255)

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

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, 5, 0.4)
    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)
    
        bodyA.update()
        bodyA.show()
        
    pygame.display.update()
    frame_rate.tick(FPS)

The Nature of Code閱讀心得筆記——使用Python實作

2024年2月29日星期四

2.6 The n-Body Problem

2.6 The n-Body Problem

Exercise 2.14

Exercise 2.15

Exercise 2.16

沒有留言:

張貼留言

4.9 Image Textures and Additive Blending

vscode

TeX

code block

2024年2月29日 星期四

2.6 The n-Body Problem

2.6 The n-Body Problem

Exercise 2.14

Exercise 2.15

Exercise 2.16

沒有留言:

張貼留言

4.9 Image Textures and Additive Blending

vscode

TeX

code block

2024年2月29日星期四