Python 繪製星球軌道教學(1)——讓一顆星球繞恆星走

in #cn7 years ago (edited)

img
太陽系示意圖(圖片來自圖庫

在開始之前,你需要的先備智識:

  • Python 基礎語法(基本I/O、迴圈等)
  • 初中數學及一點向量智識(有概念即可)

你將可以在此學到 :

  • 如何繪製單顆星球繞恆星運轉
  • 何謂數值分析
  • 微積分的概念

前言

理查‧費曼(Richard Feymann)在加州理工大學擔任教職之際,出版《費曼物理學講義》(網路上也有電子中文版,請自行搜索),書籍分述不同物理領域智識,一共三卷。希冀藉由此書出版,給予大一、大二學生清楚的觀念。

而在第一卷第九章第七節〈行星運動〉中,他曾談到行星軌道的計算:

對於吊掛在彈簧上的物體所做的震盪運動來說,以上的分析(數值分析)顯然不錯。那麼我們是否可以用類似的方法來分析行星繞太陽的運動呢?我們想知道是否可以用這個方法來得到一個近似橢圓的軌道?為了簡化問題,我們將忽略太陽本身的運動。假設一行星自某位置出發,以某一速度向前進,它沿某條曲線繞著太陽運動:我們將利用牛頓的運動定律及重力定律來分析此行星的運動,希望了解其軌跡為何。——《費曼物理學講義》

這便是此系列的主題,我們要結合物理、數學及程式來玩轉星球軌道模擬,一步步實作,探討當中的數學理論並知其物理原理,進而能模擬星球繞行恆星。

由於自身常使用 Python 撰寫程式,知道其相當好上手,不需要定義太多東西,能捨去繁瑣的語法,專注在代碼的編寫上,代碼也很容易讀;而且網路上的教學資源也蠻多,可以自己學完基礎後再來讀這教學,沒有太特別的語法與邏輯,學習曲線不會太陡,不需耗費太多心力銜接。所以在此系列,我將使用 Python 作為程式編撰的語言。

他人也可以自行用其他語言寫代碼,歡迎眾人分享自己寫的版本。而要是在學習、安裝 Python 中有問題的話,可以將問題截圖或詳細描述問題放在留言上,我會盡可能去解決你們的問題。

物理分析與推導

在此我直接講述整個推導過程,若想要更詳細的介紹就直接參閱《費曼物理學講義》的內容吧!

我們先處理二維空間上的行星繞恆星模擬,再推導到三維空間上。費曼利用一張圖(9-5,因版權問題不直接上傳)來說明星球受到的力,將力分成 Fx 和 Fy 兩道,利用相似三角形的方法,我們可以推導出下列式子:

接著使用最原始的 F = ma,進一步得到下列方程式:

依此類推,Y分量也可以推導出同等關係:

在此我們假定 GM≡1,於這種單純的一顆星球繞行恆星的情況下,只需知相對距離就好,兩者的質量不是很重要。

設定好一些參數後,我們便可以利用數值分析來模擬行星運行。但有一個參數跟數值分析息息相關,所以在編寫之前,我們還需要了解數值分析的概念才好實作。

何謂數值分析

數值分析(numerical analysis)簡而言之就是利用電腦解決一些問題。雖然不能完全搓合答案,但可以非常近似所要的解。

以這系列為例,我們會將橢圓曲線離散化,分割成有限的點,計算每個時間點星球的位置、速度、加速度。要計算這些變化,其概念就像微分一樣:

而我們的計算流程如下(記得 ϵ 是時間間隔):

時間是種連續的東西,無時無刻都在改變,時間可以不斷分割,但如此一來,資料量必然是無限。所以我們只能用一小段時間作為基本單位,將其限制在有限,改變量也隨著最小單位時間來變化。

數學定義上,ϵ 應該要多小就能有多小,也就是俗稱的無窮小。不過,電腦終究不是人類,無法處理連續的事物,我們僅能將離散盡可能化小,逼近連續。所以繪製軌道其實算是差和分而非微積分,仍有時間差需要考慮,不可能滿足無窮小的要求。要是時間差太大,繪製出的軌道會有誤差:

ϵ 越大誤差值越大,反之則越小。為了有準確的結果,我們需要一個較小的 ϵ 值,如此才能真正去模擬行星如何運行。

而費曼是這樣說的:

事實上,誤差的大小跟所採用的時間間隔 ϵ 的平方成正比。也就是說,如果我們把 ϵ 縮小一千倍,準確度會提升一百萬倍。前例中,我們採用 ϵ = 0.1 時,準確度已經達了百分之一;若要讓準確度到達十億分之一,我們還得把 ϵ 再縮小一萬倍。

但 ϵ 越小,所耗計算時間越久,不一定越小越好。有可能已經能滿足我們的條件了,卻浪費多餘時間在絲毫的差異上,所以 ϵ 只需要到達我們可容忍的範圍即可。

在這裡,我們設定 ϵ = 0.0001。作為時間間隔的值。

套件概略說明

工欲善其事,必先利其器。要運行程式前,請先確認是否安裝下列套件:

  • Numpy
  • Matplotlib

Python 本身沒有附帶這些套件,要是才剛下載下來,請安裝這些套件,詳細的安裝方法有請上網查詢。

這兩個套件常常在科學領域使用,彼此扮演不同角色。Numpy 主要處理數值運算,本身以 C 語言優化,所以效率會比普通的 Python 運算快上許多,藉此可以加速計算的速度;Matplotlib 將計算結果以圖像化呈現,有各式圖表可以使用,好比最一般的折線圖、長條圖。我們將以這個套件繪製出結果,讓人容易明白行星的運動。

程式編寫

首先我們要將需要的套件帶入,輔助接下來的模擬:

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from matplotlib import style
style.use('ggplot') #讓圖像變得好看一點,全憑個人喜好

將先前提到的 ϵ 代入,設定其為 0.0001;迭代次數則憑狀況而定,有些時候迭代次數不足,行星會走到一半就結束,原本要的橢圓就會破一個洞,不怎麼完美。一般情形下,我們會將星球的初始位置設定在 1.5 以內,如此才不會運行太久。事實上,我們模擬的目的主要是看其軌道的型態為何,彼此間的距離其實沒什麼關係,在相等比例下,繪製出的軌道基本上會相同,沒什麼必要讓它繞很大一圈。

以下是我們設定的兩個超參數:

epsilon = 0.0001 #時間間隔
n = 100000 #迭代次數

思考一下,星球本身含有那些屬性需要我們考慮呢?位置、速度、加速度,這些數值我們需要記錄下來,為運算必備之元素:

position = np.zeros((n, 2)) #紀錄每個時間點的距離
velocity = np.zeros((2)) #紀錄當下的速度
acceleration = np.zeros((2)) #紀錄當下的加速度

除卻要記錄這些數據外,我們還需要初始值,如此才能知道該如何運作。那模擬星球繞行恆星需要輸入幾個參數呢?位置是一個,而速度則是另外一個;此外,我們在二維空間上運行,所以位置與速度皆有 X、Y 軸方向之分,各兩個參數,於是一共要輸入4筆資料。

輸入資料:
  • X初始座標值
  • Y初始座標值
  • X分量初始速度
  • Y分量初始速度

有關輸入參數的代碼如下:

position[0, 0] = float(input("輸入X初始座標值: ")) 
position[0, 1] = float(input("輸入Y初始座標值: "))
velocity[0] = float(input("輸入X分量初始速度: "))
velocity[1] = float(input("輸入Y分量初始速度: "))

得到參數後,我們便要開始一連串的計算。還記得上頭所寫的方程式嗎?現在我們要將它們轉化為程式碼,讓電腦替我們計算。我們依照費曼書中提到的方式,先算加速度與速度。流程大抵如下:

  1. 需要知道星球距離恆星多遠:求 r。
  2. 需要計算兩者間距離三次方的倒數: 求 1 / r3
  3. 需要計算星球的速度:求 v。

代碼長得像這樣子:

r = norm(position[0]) #計算離太陽距離
r3 = 1 / (r**3) #計算 1/r^3
acceleration = position[0] * -r3 #計算星球加速度
velocity += acceleration * epsilon * 0.5 #計算星球速度

接下來,我們可以開始一連串的運算了。每回合流程都一樣,我們可以用個迴圈來處理。這是最核心的代碼,等一下移植到三維空間,也是使用下列的代碼,原封不動:

i = 1
while(i < n): #重複n次
    position[i] = position[i-1] + velocity * epsilon #計算更新後的星球位置
    r = norm(position[i]) #計算恆星與星球間的距離
    r3 = 1 / (r**3) #計算 1/r^3
    acceleration = position[i] * -r3 #計算星球加速度
    velocity += acceleration * epsilon #計算星球速度
    i += 1

算完後,我們要繪製出計算後的結果:

plt.plot(position[:, 0], position[:, 1], 'b-')
plt.plot(0, 0, "ro")
plt.show()

整個星球運動模擬程式就大功告成,大致上並不複雜,三、四十行代碼就能搞定。

讓我們把《費曼物理學講義》上的參數(vx=0,vy=1.63,x=0.5,y=0)代入測試,看看會得出什麼樣的結果來:

看來還不錯,行星繞著恆星的軌道路徑是一顆豐滿的橢圓,就像蛋一樣。

拓展到三維

其實並不難把行星路徑模擬搬到三維空間運行,所有原理一樣;我們依舊能在各個分量找到它們的相似三角形,再用相同推導方式求出方程式來。唯一不同之處是計算距離的公式需要修改:

但我們使用的 Numpy 套件,其函數本身就可以計算高維度空間的向量長度,也就是求出它們的模長,所以大致上不需做任何修改。需要更動的地方是 Matplotlib 繪圖,原先它只能繪製 2D 圖形,現在我們要把它弄成可以繪製 3D 圖形,整體代碼如下:

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import style
style.use('ggplot')
fig = plt.figure()
Ax = fig.gca(projection='3d')
epsilon = 0.0001 #間隔
n = 100000 #迭代次數
position = np.zeros((n, 3))
velocity = np.zeros((3))
acceleration = np.zeros((3))
position[0, 0] = float(input("輸入X初始座標值: ")) 
position[0, 1] = float(input("輸入Y初始座標值: "))
position[0, 2] = float(input("輸入Z初始座標值: "))
velocity[0] = float(input("輸入X分量初始速度: "))
velocity[1] = float(input("輸入Y分量初始速度: "))
velocity[2] = float(input("輸入Z分量初始速度: "))
r = norm(position[0])
r3 = 1 / (r**3) 
acceleration = position[0] * -r3
velocity += acceleration * epsilon * 0.5
i = 1
while(i < n): #核心代碼不變
    position[i] = position[i-1] + velocity * epsilon
    r = norm(position[i])
    r3 = 1 / (r**3)
    acceleration = position[i] * -r3
    velocity += acceleration * epsilon
    i += 1
Ax.plot(position[:, 0], position[:, 1], position[:, 2], 'b-')
Ax.scatter(0, 0, 0, c="r")
for angle in range(0, 360): #360度旋轉一圈
    Ax.view_init(30, angle)
    plt.draw()
    plt.pause(.001)

下面是我隨意代參數進去後所得的結果,大家可以自行試試:

有關星球軌道的基礎繪製就在此暫時告一段落,但旅程尚未終了,還有許多進階的事物等著我們探索。

不過,在此之前,你們可以試著修改代碼使其臻於完善,或者就純粹把玩它看看能玩出甚麼新花樣。總之,我期待眾人能有所啟發,能有一些新想法,精益求精,開始創造自己的小宇宙吧!

思考與改進

  • 會不會覺得行星已經繞一圈後仍繼續計算有點多餘?有沒有方法可以解決這個問題呢?
  • 行星在一定條件下,會擺脫恆星的重力,逕自向外飛去而不復返。是否只能等到計算完才能結束呢?抑或是有特殊的算法可以知道行星是否會脫離?
小弟才學疏淺,內容難免有所紕漏,盼望眾人可以容忍並提醒我更正。
Sort:  
Loading...

这个覆盖面太广了!!!根本没法吸收知识啊 好厉害的感觉

這的確覆蓋很多東西。我原意是想要將各式學科結合在一起,而不要死板板地限制在一個領域之中。在學校,老師大都侷限於自己專門的領域中授課,各個科目沒什麼交集。這樣感覺不怎麼好,就好比一隻魚被大卸八塊卻不鮮美可口,成了屍塊一般。老實說,學科並非涇渭分明而有個明確的界線,若肯努力去探索,必定可以找到它們的異同處。

台灣最近在鼓吹程式教育,希望整個教育體系能好好納入程式設計。這想法本身不錯,自己也樂觀如是發展,但擔心會流於形式,而失去其中的真諦。

所以我才寫這系列,想作為一個範本供他們參考。這一篇文章本身就是高中領域內的東西來做發揮,行星運行、牛頓力學、向量,沒有超出範圍,如果是好好學習的高中生,應並不難理解當中的理論。這樣一來,在授課的同時便可以結合程式設計,一同學習,不一定要特別開一門資訊課來教導。

不過,我是期望大家都可以學習,好好地體會當中的奧妙啦!只可惜這系列難度有點高,算是一個遺憾。

Great post and very informative. @etzel Have learnt a lot from it.

我對於這個就真的完全沒有概念 @@
不過很歡迎你可以來這裡分享你的所知喔! :)

沒關係,總會有人懂。

我確實很想在這裡分享我所學,有好多東西可以說喔!

可以喔!! 這裡有很多不同專業領域的朋友! :)

屌,很高興又在這平台遇到臺灣人!

應該不難遇到吧!呵!

最近有一對退休夫妻也來這裡了,同樣是台灣人,不知你有沒有看到。

其實也沒這麼容易喔
很高興在這裡認識你,歡迎你加入。 :D

他們的帳號是什麼呢?

@tungjungchen@helenliu兩位,女兒邀他們來這。

有看到他們,很有意思的一家人