用Python搓一個太陽系

日地月三體

所謂三體，就是三個物體在重力作用下的運動，由于三點共面，所以三個質點僅在重力作用下的運動軌跡也必然無法逃離平面，

三體運動所遵循的規律就是古老而經典的萬有引力

F ? = G m i m j r 2 e ? r \vec F=\frac{Gm_im_j}{r^2}\vec e_r F =r2Gmi?mj??e r?

則對于 m i m_i mi?而言，

m i d v ? i d t = G m i m j r i j 3 r ? i j m_i\frac{\text d\vec v_i}{\text dt}=\frac{Gm_im_j}{r_{ij}^3}\vec r_{ij} mi?dtdv i??=rij3?Gmi?mj??r ij?

且

d r ? i d t = v ? i \frac{\text d\vec r_i}{\text dt}=\vec v_i dtdr i??=v i?

將其寫為差分形式

v ? i = ∑ j =? i G m j r i j 3 r ? i j d t r ? i = v ? i d t \begin{aligned} \vec v_i&=\sum_{j\not=i}\frac{Gm_j}{r_{ij}^3}\vec r_{ij}\text dt\\ \vec r_i&= \vec v_i\text dt \end{aligned} v i?r i??=j?=i∑?rij3?Gmj??r ij?dt=v i?dt?

由于我們希望觀察三體運動的復雜形式，而不關系其隨對應的宇宙星體，所以不必考慮單位制，將其在二維平面坐標系中拆分，令 v ? = ( u , v ) \vec v=(u,v) v =(u,v)，則

u i + = ∑ j =? i G m j ( x j ? x i ) d t ( x i ? x j ) 2 + ( y i ? y j ) 2 3 v i + = ∑ j =? i G m j ( y j ? y i ) d t ( x i ? x j ) 2 + ( y i ? y j ) 2 3 x i + = u ? i d t y i + = v ? i d t \begin{aligned} u_i&+=\sum_{j\not=i}\frac{Gm_j(x_j-x_i)\text dt}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\\ v_i&+=\sum_{j\not=i}\frac{Gm_j(y_j-y_i)\text dt}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\\ x_i&+= \vec u_i\text dt\\ y_i&+= \vec v_i\text dt \end{aligned} ui?vi?xi?yi??+=j?=i∑?(xi??xj?)2+(yi??yj?)2 ?3Gmj?(xj??xi?)dt?+=j?=i∑?(xi??xj?)2+(yi??yj?)2 ?3Gmj?(yj??yi?)dt?+=u i?dt+=v i?dt?

太陽、地球和月亮就是一個典型的三體系統，其中太陽質量為 1.989 × 1 0 30 k g 1.989×10^{30}kg 1.989×1030kg，地球質量為 5.965 × 1 0 24 k g 5.965×10^{24}kg 5.965×1024kg，月球質量為 7.342 ? 1 0 22 k g 7.342?10^{22}kg 7.342?1022kg，萬有引力常數為 G = 6.67 × 1 0 ? 11 N ? m 2 / k g 2 G=6.67×10^{-11}N·m2/kg^2 G=6.67×10?11N?m2/kg2，地月距離為 3.8 × 1 0 8 m 3.8\times10^8m 3.8×108m；日地距離為 1.5 × 1 0 11 m 1.5\times10^{11}m 1.5×1011m；地球公轉速度為 28.8 k m / s 28.8km/s 28.8km/s；月球公轉速度為 1 k m / s 1km/s 1km/s，則各引數初始化為

#后續代碼主要更改這里的引數
m = [1.33e20,3.98e14,4.9e12]
x = np.array([0,1.5e11,1.5e11+3.8e8])
y = np.array([0,0,0])
u = np.array([0,0,0])
v = np.array([0,2.88e4,1.02e3])

由于地月之間的距離相對于日地距離太近，所以在畫圖的時候將其擴大100倍，得到影像

盡管存在誤差，但最起碼看到了地球圍繞太陽轉，月球圍繞地球轉，，，代碼為

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

m = [1.33e20,3.98e14,4.9e12]
x = np.array([0,1.5e11,1.5e11+3.8e8])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.88e4+1.02e3])

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(xlim=(-2e11,2e11),ylim=(-2e11,2e11))
ax.grid()

trace0, = ax.plot([],[],'-', lw=0.5)
trace1, = ax.plot([],[],'-', lw=0.5)
trace2, = ax.plot([],[],'-', lw=0.5)

pt0, = ax.plot([x[0]],[y[0]] ,marker='o')
pt1, = ax.plot([x[0]],[y[0]] ,marker='o')
pt2, = ax.plot([x[0]],[y[0]] ,marker='o')

k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'

N = 1000
dt = 36000
ts =  np.arange(0,N*dt,dt)/3600/24
xs,ys = [],[]
for _ in ts:
    x_ij = (x-x.reshape(3,1))
    y_ij = (y-y.reshape(3,1))
    r_ij = np.sqrt(x_ij**2+y_ij**2)
    for i in range(3):
        for j in range(3):
            if i!=j :
                u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)
                v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)
    x += u*dt
    y += v*dt
    xs.append(x.tolist())
    ys.append(y.tolist())

xs = np.array(xs)
ys = np.array(ys)

def animate(n):
    trace0.set_data(xs[:n,0],ys[:n,0])
    trace1.set_data(xs[:n,1],ys[:n,1])
    #繪圖時的地月距離擴大100倍，否則看不清
    tempX2S = xs[:n,1]+100*(xs[:n,2]-xs[:n,1])
    tempY2S = ys[:n,1]+100*(ys[:n,2]-ys[:n,1])
    trace2.set_data(tempX2S,tempY2S)

    pt0.set_data([xs[n,0]],[ys[n,0]])
    pt1.set_data([xs[n,1]],[ys[n,1]])
    tempX = xs[n,1]+100*(xs[n,2]-xs[n,1])
    tempY = ys[n,1]+100*(ys[n,2]-ys[n,1])
    pt2.set_data([tempX],[tempY])

    k_text.set_text(textTemplate % ts[n])
    return trace0, trace1, trace2, pt0, pt1, pt2, k_text

ani = animation.FuncAnimation(fig, animate, 
    range(N), interval=10, blit=True)

plt.show()
ani.save("3.gif")

日地火

m = [1.33e20,3.98e14,4.28e13]
x = np.array([0,1.5e11,2.28e11])
y = np.array([0.0,0,0])
u = np.array([0.0,0,0])
v = np.array([0,2.88e4,2.4e4])

### 由于火星離地球很遠，所以不必再改變尺度
def animate(n):
    trace0.set_data(xs[:n,0],ys[:n,0])
    trace1.set_data(xs[:n,1],ys[:n,1])
    trace2.set_data(xs[:n,2],ys[:n,2])

    pt0.set_data([xs[n,0]],[ys[n,0]])
    pt1.set_data([xs[n,1]],[ys[n,1]])
    pt2.set_data([xs[n,2]],[ys[n,2]])

    k_text.set_text(textTemplate % ts[n])
    return trace0, trace1, trace2, pt0, pt1, pt2, k_text

得到

這個運動要比月球的運動簡單得多——前提是開上帝視角，俯瞰太陽系，如果站在地球上觀測火星的運動，那么這個運動可能相當帶感

所以這都能找到規律，托勒密那幫人也真夠有才的，

太陽系

由于太陽和其他星體之間的質量相差懸殊，所以太陽系內的多體運動，都將退化為二體問題，甚至如果把太陽當作不動點，那就成了單體問題了，

盡管如此，我們還是盡可能地模仿一下太陽系的運動情況

除了水星偏心率為0.2，對黃道面傾斜為7°之外，其余行星的偏心率皆小于0.1，且對黃道面傾斜普遍小于4°，由于水星的軌道太小，偏不偏心其實都不太看得出來，所以就當它是正圓也無所謂了，最后得圖

au,G,RE,ME = 1.48e11,6.67e-11,1.48e11,5.965e24

m = np.array([3.32e5,0.055,0.815,1,
              0.107,317.8,95.16,14.54,17.14])*ME*6.67e-11
r = np.array([0,0.387,0.723,1,1.524,5.203,
              9.537,19.19,30.7])*RE
theta = np.random.rand(9)*np.pi*2

x = r*np.cos(theta)
y = r*np.sin(theta)

v = np.array([0,47.89,35.03,29.79,
              24.13,13.06,9.64,6.81,5.43])*1000
u = -v*np.sin(theta)
v = v*np.cos(theta)

name = "solar.gif"

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(xlim=(-31*RE,31*RE),ylim=(-31*RE,31*RE))
ax.grid()

traces = [ax.plot([],[],'-', lw=0.5)[0] for _ in range(9)]
pts = [ax.plot([],[],marker='o')[0] for _ in range(9)]

k_text = ax.text(0.05,0.85,'',transform=ax.transAxes)
textTemplate = 't = %.3f days\n'

N = 500
dt = 3600*50
ts =  np.arange(0,N*dt,dt)
xs,ys = [],[]
for _ in ts:
    x_ij = (x-x.reshape(len(m),1))
    y_ij = (y-y.reshape(len(m),1))
    r_ij = np.sqrt(x_ij**2+y_ij**2)
    for i in range(len(m)):
        for j in range(len(m)):
            if i!=j :
                u[i] += (m[j]*x_ij[i,j]*dt/r_ij[i,j]**3)
                v[i] += (m[j]*y_ij[i,j]*dt/r_ij[i,j]**3)
    x += u*dt
    y += v*dt
    xs.append(x.tolist())
    ys.append(y.tolist())

xs = np.array(xs)
ys = np.array(ys)

def animate(n):
    for i in range(9):
        traces[i].set_data(xs[:n,i],ys[:n,i])
        pts[i].set_data(xs[n,i],ys[n,i])
    k_text.set_text(textTemplate % (ts[n]/3600/24))
    return traces+pts+[k_text]

ani = animation.FuncAnimation(fig, animate, 
    range(N), interval=10, blit=True)

plt.show()
ani.save(name)

由于外圈的行星軌道又長速度又慢，而內層的剛好相反，所以這個圖很難兼顧，觀感上也不太好看，

如果只畫出木星之前的星體，順便加上小行星帶，可能會好一些，

通過這個圖就能看出來，有一顆小行星被木星彈了過來，直沖沖地向地球趕來，幸好又被太陽彈了出去，可見小行星還是挺危險的，好在這只是個假想圖，