行星運動之數值模擬

以Euler method計算行星運動

行星運動的動力就是牛頓的重力定律，而運動的法則就是牛頓第二運動定律。但是各位同學要特別留意到，這一個運動是一個平面的運動，也就是在二維的運動，與拋體運動相似，都是在平面上的運動。因此我們在做差分近似的時候，需要使用兩個座標變數(x,y)和兩個速度分量(vx,vy)，當然也有兩個加速度分量(ax,ay)。而對應的牛頓運動方程式也是一個向量方程式：\[\vec{F}=m\vec{a}, \,\,\, \vec{a}=\frac{\vec{F}}{m}\] \[\vec{a}=\frac{d\vec{v}}{dt}=\frac{\vec{F}}{m}\] 而行星運動時所遭遇到的作用力就是行星與太陽之間的重力，根據牛頓的重力定律： \[\vec{F}=-\frac{GM_s m}{r^3}\vec{r}, \] 其中\(M_s\)是太陽的質量，\(m\)是行星的質量，\(r\)太陽與行星之間的距離，\(G=6.67 \times 10^{-11}\)是重力的宇宙常數。因此行星在太陽的重力場中運動的方程式如下： \[ \frac{d\vec{r}}{dt}=\vec{v} \] \[ \frac{d\vec{v}}{dt}=\vec{a}=\frac{\vec{F}}{m}=-\frac{GM_s}{r^3}\vec{r} \]
運動的差分近似(Euler method) \[ \vec{r}(t+dt) \simeq \vec{r}(t)+\vec{v}(t)dt \] \[ \vec{v}(t+dt) \simeq \vec{v}(t)+\vec{a}(t)dt \] 為了簡化我們整個物理問題的公式形式，我們要選擇一個恰當的單位系統，如果將軌道運動視為一個圓週運動，那麼我們可以根據圓週運動的向心加速度的公式和週期的公式得到下面的結果： \[\frac{4\pi^2 m r}{T^2}=\frac{GM_s m}{r^2},\] \[\frac{r^3}{T^2}=\frac{GM_s}{4\pi^2} \,\,\,\, \text{(Kepler's 3rd law)}\] 這就是克卜勒的行星運動第三定律，行星運動的週期(T)的平方與行星運動的平均距離(r)的3次方成正比。將這個公式運用在我們的地球公轉上我們將會得到，\(T=1年，r=1AU\)，AU=天文單位=地球與太陽之間的距離=\(1.5 \times 10^{11}\) 公尺。如果我們將時間的單位訂為年，距離的單位定為天文單位，那麼這個公式將簡化為 \[\frac{1^3}{1^2}=\frac{GM_s}{4\pi^2}, \,\,\, GM_s=4\pi^2\] 因此在這個單位系統中\(GM_s=4\pi^2\)，地球在某一個時刻的坐標為\(\vec{r}=(1,0)\)，而這個時刻的速度為\(\vec{v}=(0,2\pi) \,\, v=\frac{2\pi r}{T}=\frac{2\pi \cdot 1}{1}=2\pi\)。



---

程式樣本1,1個行星運動(Euler method)

import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("Agg")
import numpy as np

m=1.0; GM=4.*np.pi**2
R0=1.; V0=np.pi*2.
r1=np.array([R0,0.])
v1=np.array([0., V0])
r1n=np.linalg.norm(r1)
a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
xp1=[r1[0]]; yp1=[r1[1]]
t=0.; dt=0.01
tp=[t]
while t < 2.:
    print ('%11.5f'*6 %(t,r1[0],r1[1],v1[0],v1[1],r1n))
    #print '%9.4f'*7 %(t,dr,drx,r1[0],r1[1],r2[0],r2[1])
    v1n=np.linalg.norm(v1)
    r1n=np.linalg.norm(r1)
    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    r1+=v1*dt  #Euler
    v1+=a1*dt
    #r1+=v1*dt #Euler-Cromer
    t+=dt
    tp.append(t)
    xp1.append(r1[0])
    yp1.append(r1[1])
plt.plot(xp1,yp1)
plt.axis('equal')
plt.savefig("orbital-1.png")
plt.close()







---

程式樣本2,2個行星運動(Euler method)

# coding=Big5
import numpy as np
import matplotlib.pyplot as plt
m=1.0; GM=4.*np.pi**2
R0=1.; V0=np.pi*2.

r1=np.array([R0,0.])
v1=np.array([0., V0])
r1n=np.linalg.norm(r1)    
a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
xp1=[r1[0]]; yp1=[r1[1]]

r2=np.array([R0,0.])
v2=np.array([0., V0*0.8])
r2n=np.linalg.norm(r2)    
a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
xp2=[r2[0]]; yp2=[r2[1]]

t=0.; dt=0.01
tp=[t]
while t < 2.:
    print '%10.5f'*7 %(t,r1[0],r1[1],r2[0],r2[1],r1n,r2n)
    #print '%9.4f'*7 %(t,dr,drx,r1[0],r1[1],r2[0],r2[1])

    v1n=np.linalg.norm(v1)
    r1n=np.linalg.norm(r1)
    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    r1+=v1*dt
    v1+=a1*dt
    #r1+=v1*dt

    v2n=np.linalg.norm(v2)
    r2n=np.linalg.norm(r2)
    a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
    r2+=v2*dt
    v2+=a2*dt
    #r2+=v2*dt
    
    t+=dt
    tp.append(t)
    xp1.append(r1[0]); yp1.append(r1[1])    
    xp2.append(r2[0]); yp2.append(r2[1])
    
plt.plot(xp1,yp1)
plt.axis('equal')
plt.plot(xp2,yp2)
plt.show()







Euler-Cromer method

當我們在使用Euler method的時候我們用了兩次差分近似，一次是位置的差分近似，另外一個是速度的差分近似。當我們在使用他的時候，原則上應該都是用前一個時刻的位置與速度，因此在程式書寫的時候，到底哪一個差分近似先執行都沒有影響。但是如果我們先將速度更新，再用更新後的速度來做位置的更新，那麼這個方法就稱為Euler-Cromer method. 也就是說如果我們要使用Euler

    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    r1+=v1*dt
    v1+=a1*dt

那麼這樣的順序是對的，因為r1的更新時所使用的速度v1，就是時間為t的時候的速度；v1更新時所使用的加速度，也是在時間為t的時候的加速度a1。Euler-Cromer method是將位置的更新與速度的更新的兩條指令對調：

    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    v1+=a1*dt
    r1+=v1*dt

如此一來雖然速度v1的更新，仍然是使用前一個時刻的加速度a1，但是位置r1的更新，所使用的速度並前一時刻的速度，而是已經更新過後的速度。所以如果我們的數值計算是用這種方式來進行的話，那麼將是兩種更新的混合。原則上我們並不知道哪一個方式會得到比較精確的結果，但是在實際應用中，似乎Euler-Cromer method都具有較高的精確度，例如下面的樣本程式就是比較Euler和Euler-Cromer計算一個行星運動的程式，我們從這個程式執行得到的數據作比較，可以明顯發現Euler-Cromer確實有相當大幅度的改進。



---

程式樣本3,行星運動(Euler-Cromer method)

import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("Agg")
import numpy as np

m=1.0; GM=4.*np.pi**2
R0=1.; V0=np.pi*2.

r1=np.array([R0,0.])
v1=np.array([0., V0])
r1n=np.linalg.norm(r1)
a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
xp1=[r1[0]]; yp1=[r1[1]]

r2=np.array([R0,0.])
v2=np.array([0., V0*1.])
r2n=np.linalg.norm(r2)
a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
xp2=[r2[0]]; yp2=[r2[1]]

t=0.; dt=0.01
tp=[t]
while t < 2.:
    print ('%10.5f'*7 %(t,r1[0],r1[1],r2[0],r2[1],r1n,r2n))

    v1n=np.linalg.norm(v1)
    r1n=np.linalg.norm(r1)
    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    r1+=v1*dt
    v1+=a1*dt
    #r1+=v1*dt

    v2n=np.linalg.norm(v2)
    r2n=np.linalg.norm(r2)
    a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
    #r2+=v2*dt
    v2+=a2*dt
    r2+=v2*dt

    t+=dt
    tp.append(t)
    xp1.append(r1[0]); yp1.append(r1[1])
    xp2.append(r2[0]); yp2.append(r2[1])

plt.plot(xp1,yp1,label=r"Euler")
plt.plot(xp2,yp2,label=r"Euler-Cromer")
plt.axis('equal')
plt.legend()
plt.savefig("orbital-3.png")
plt.close()







Runge-Kutta 2nd order method

我們在前兩個單元已經可以看到RK2的方法在改進Euler method的精確度上面，確實有非常大的進展，現在我們也在軌道運動問題當中引入RK2的數值方法。
• 對單一一個微分方程\[\frac{dx}{dt}=f(x,t)\]
RK2的做法如下：\[x(t+\tau)=x(t)+\tau f^*(x^*,t^*)\]
• 在自由落體運動中
我們已經使用過RK2的方法，其中面臨有兩個微分方程式需要處理，一個是位置(\(\frac{dy}{dt}=v\))另外一個是速度(\(\frac{dv}{dt}=a\))，所使用的公式如下：
根據RK2的原理，首先要算出在\(t^*=t+dt/2\)的\[y^*=y(t)+v(t)*dt/2,\] \[ v^*=v(t)+a(t)*dt/2,\] \[ a^*=a(t^*,v^*)=g-b/m*v^*\]，再利用這些中介的數據求取新的位置與速度： \[y(t+dt)+=v^**dt,\] \[ v(t+dt)+=a^**dt \] \[t^*=t+\tau/2\]
• 對於現在我們要處理的行星軌道運動問題，
基本上也是有兩個微分方程，位置與速度兩個微分方程，但是更複雜的是行星在平面上運動，每一個方程式都有兩個分量(x,y,vx,vy)，因此RK2的公式應為： \[ \vec{r}(t+\tau)=\vec{r}(t) + \tau \vec{v}^*\] \[ \vec{v}(t+\tau)=\vec{v}(t) + \tau \vec{a}^*\] 其中的中介速度與加速度向量為 \[ \vec{r}^*=\vec{r}(t) + \frac{\tau}{2}\vec{v}(t) \] \[ \vec{v}^*=\vec{v}(t) + \frac{\tau}{2}\vec{a}(t) \] \[ \vec{a}^*=\vec{a}(t^*, \vec{r}^*, \vec{v}^*)=-\frac{GM_s}{{r^*}^3}\vec{r}^* \] \[t^*=t+\tau/2\]
根據這些描述我們在下面提供的一個樣本程式，讓同學們可以參考如何使用RK2的方法來實現行星運動的數值模擬，並且你也可以透過數據間的比較理解到RK2確實在精確度上提供了很好的改善。



---

程式樣本4,行星運動(RK2 method)

import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("Agg")
import numpy as np

m=1.0; GM=4.*np.pi**2
R0=1.; V0=np.pi*2.

r1=np.array([R0,0.])
v1=np.array([0., V0])
r1n=np.linalg.norm(r1)
a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
xp1=[r1[0]]; yp1=[r1[1]]

r2=np.array([R0,0.])
v2=np.array([0., V0*1.])
r2n=np.linalg.norm(r2)
a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
xp2=[r2[0]]; yp2=[r2[1]]

t=0.; dt=0.01; dt2=dt/2.
tp=[t]
while t < 2.:
    print ('%10.5f'*7 %(t,r1[0],r1[1],r2[0],r2[1],r1n,r2n))

    v1n=np.linalg.norm(v1)
    r1n=np.linalg.norm(r1)
    a1=np.array([-GM/r1n**3*r1[0],-GM/r1n**3*r1[1]])
    r1+=v1*dt
    v1+=a1*dt
    #r1+=v1*dt

    v2n=np.linalg.norm(v2)
    r2n=np.linalg.norm(r2)
    a2=np.array([-GM/r2n**3*r2[0],-GM/r2n**3*r2[1]])
    ts=+dt2
    r2s=r2+v2*dt2
    v2s=v2+a2*dt2
    r2sn=np.linalg.norm(r2s)
    a2s=np.array([-GM/r2sn**3*r2s[0],-GM/r2sn**3*r2s[1]])

    r2+=v2s*dt
    v2+=a2s*dt
    #r2+=v2*dt

    t+=dt
    tp.append(t)
    xp1.append(r1[0]); yp1.append(r1[1])
    xp2.append(r2[0]); yp2.append(r2[1])
plt.plot(xp1,yp1,label=r"Euler")
plt.plot(xp2,yp2,label=r"RK2")
plt.axis('equal')
plt.legend()
plt.savefig("orbital-4.png")
plt.close()