【Python】有三颗恒星的三体人很难不产生超能力

三体问题

经同学提出，我才意识到，原来三体人有三颗恒星……也就意味着可能三体星人连个稳定的恒星轨道都没有，悲惨指数直线上升。

但就拉格朗日方程而言，却并不困难。

设 m i , i ∈ { 0 , 1 , 2 , 3 } m_i, i\in\{0,1,2,3\} mi​,i∈{0,1,2,3}表示四颗星体，则对任意星体 i i i而言，其动能为

T i = 1 2 m i ( x ˙ 2 + y ˙ 2 ) T_i=\frac1 2m_i(\dot x^2+\dot y^2) Ti​=21​mi​(x˙2+y˙​2)

势能为

V i = − ∑ j ≠ i G m i m j ( x i − x j ) 2 + ( y i − y j ) 2 V_i=-\sum_{j\not=i}\frac{Gm_im_j}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}} Vi​=−j​=i∑​(xi​−xj​)2+(yi​−yj​)2 ​Gmi​mj​​

拉格朗日量为 L = T − V L=T-V L=T−V，根据拉格朗日方程

d d t ∂ L ∂ r ˙ − ∂ L ∂ r = 0 ， r = x i , y i \frac{\text d}{\text dt}\frac{\partial L}{\partial\dot r}-\frac{\partial L}{\partial r}=0，r=x_i,y_i dtd​∂r˙∂L​−∂r∂L​=0，r=xi​,yi​

有

x ¨ i = − ∑ i ≠ j G m j ( x i − x j ) ( x i − x j ) 2 + ( y i − y j ) 2 3 y ¨ i = − ∑ i ≠ j G m j ( y i − y j ) ( x i − x j ) 2 + ( y i − y j ) 2 3 \ddot x_i=-\sum_{i\not =j}\frac{Gm_j(x_i-x_j)}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\ \ddot y_i=-\sum_{i\not =j}\frac{Gm_j(y_i-y_j)}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3} x¨i​=−i​=j∑​(xi​−xj​)2+(yi​−yj​)2 ​3Gmj​(xi​−xj​)​y¨​i​=−i​=j∑​(xi​−xj​)2+(yi​−yj​)2 ​3Gmj​(yi​−yj​)​

则令state以如下顺序编排， s = x 0 , x ˙ 0 , y 0 , y ˙ 0 , x 1 , x ˙ 1 … s=x_0,\dot x_0, y_0, \dot y_0, x_1,\dot x_1\dots s=x0​,x˙0​,y0​,y˙​0​,x1​,x˙1​…，则 s 4 i = x i , s 4 i + 1 = x ˙ i , s 4 i + 2 = y i , s 4 i + 3 = y ˙ i s_{4i}=x_i, s_{4i+1}=\dot x_i, s_{4i+2}=y_i, s_{4i+3}=\dot y_i s4i​=xi​,s4i+1​=x˙i​,s4i+2​=yi​,s4i+3​=y˙​i​。

则列出微分方程如下

def derivs(state, t):
    N = int(len(state)/4)
    dydx = np.zeros_like(state)
    for i in range(N*2):
        dydx[i*2] = state[i*2+1]
    for i in range(N):
        dydx[i*4+1] = 0
        dydx[i*4+3] = 0
        for j in range(N):
            if i==j:continue
            dx = state[i*4]-state[j*4]
            dy = state[i*4+2]-state[j*4+2]
            L = np.sqrt(dx**2+dy**2)**3
            dydx[i*4+1] -= G * m[j] * dx / L
            dydx[i*4+3] -= G * m[j] * dy / L
    return dydx

由于三体运动过于放荡不羁，故而随机生成的三体几乎很快就分道扬镳了，所以接下来选择适当位置和重量的三颗恒星。

且令万有引力常数以年为时间单位

G = 4.98 × 1 0 − 10 k m 3 d − 2 k g − 1 G=4.98\times10^{-10} km^3d^{-2}kg^{-1} G=4.98×10−10km3d−2kg−1

令恒星质量在 1 0 30 k g 10^{30}kg 1030kg的量级，空间距离在 1 0 11 k m 10^{11}km 1011km量级。

行星质量在 1 0 25 10^{25} 1025量级。

由于质量相差过大，所以假定行星质量为0也是可以的。

为了让恒星三体尽量稳定，在生成质量和初始坐标之后，令其初速度约等于稳定三体运动的速度。

首先，星体质量为 m i m_i mi​，坐标为 ( X i , Y i ) (X_i,Y_i) (Xi​,Yi​)，则其重心坐标为

x g = ∑ i m i X i ∑ m i ， y g = ∑ i m i Y i ∑ m i x_g = \frac{\sum_im_iX_i}{\sum m_i}，y_g = \frac{\sum_im_iY_i}{\sum m_i} xg​=∑mi​∑i​mi​Xi​​，yg​=∑mi​∑i​mi​Yi​​

如将坐标系移动到 ( x g , y g ) (x_g,y_g) (xg​,yg​)，则新坐标系下 x i = X i − x g , y i = Y i − x g x_i=X_i-x_g, y_i=Y_i-x_g xi​=Xi​−xg​,yi​=Yi​−xg​。

则对 m i m_i mi​而言，其运动的半径与加速度分别为为

r i = x i 2 + y i 2 x ¨ i = ∑ j ≠ i G m j ( x j − x i ) ( x i − x j ) 2 + ( y i − y j ) 2 3 y ¨ i = ∑ j ≠ i G m j ( y j − y i ) ( x i − x j ) 2 + ( y i − y j ) 2 3 ω i = x ¨ i 2 + y ¨ i 2 r i \begin{aligned} r_i&=\sqrt{x_i^2+y_i^2}\ \ddot x_i&=\sum_{j\not=i}\frac{Gm_j(x_j-x_i)}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\ \ddot y_i&=\sum_{j\not=i}\frac{Gm_j(y_j-y_i)}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2}^3}\ \omega_i&=\sqrt{\frac{\sqrt{\ddot x_i^2+\ddot y_i^2}}{r_i}} \end{aligned} ri​x¨i​y¨​i​ωi​​=xi2​+yi2​ ​=j​=i∑​(xi​−xj​)2+(yi​−yj​)2 ​3Gmj​(xj​−xi​)​=j​=i∑​(xi​−xj​)2+(yi​−yj​)2 ​3Gmj​(yj​−yi​)​=ri​x¨i2​+y¨​i2​ ​​ ​​

如果 ω i \omega_i ωi​不相等，那么每个星体的角速度不同，则运行之后会马上打破现有的局面，从而进入不稳定的状态。

则其初始角度和速度为

cos ⁡ θ = x i r i sin ⁡ θ = y i r i u = x ˙ = − ω i r i sin ⁡ θ = − ω i y i v = y ˙ = ω i r i cos ⁡ θ = ω i x i \begin{aligned} \cos\theta&=\frac{x_i}{r_i}\ \sin\theta&=\frac{y_i}{r_i}\ u=\dot x&=-\omega_ir_i\sin\theta&=&-\omega_iy_i\ v=\dot y&=\omega_ir_i\cos\theta&=&\omega_ix_i \end{aligned} cosθsinθu=x˙v=y˙​​=ri​xi​​=ri​yi​​=−ωi​ri​sinθ=ωi​ri​cosθ​==​−ωi​yi​ωi​xi​​

得到下图，其中轨迹比较细的那个是行星……

则其初始化方法为

# 用于初始化星体的质量和位置
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from scipy import integrate

np.random.seed(42)
G = 4.98e-10
m = np.random.rand(4)*10e30
m[0] /= 1e5     #行星质量
x0 = np.random.rand(4)*1e9
x0[0] *= 2      #让行星尽量离他们三颗恒星远一点
y0 = np.random.rand(4)*1e9
y0[0] *= 2
M = np.sum(m)       #总质量
# 计算质心，并以质心为原点
x0 -= np.sum(m*x0)/M
y0 -= np.sum(m*y0)/M
r = np.sqrt(x0**2+y0**2)
a = np.zeros(4)
b = np.zeros(4)
for i in range(4):
    for j in range(4):
        if i==j : continue
        dx = x0[i]-x0[j]
        dy = y0[i]-y0[j]
        L = np.sqrt(dx**2+dy**2)**3
        gm = G * m[j]
        a[i] += gm * dx / L
        b[i] += gm * dy / L

om = np.sqrt(np.sqrt(a**2+b**2)/r)
u0 = -om*y0
v0 = om*x0

绘图代码为

state = np.zeros(16)
for i in range(4):
    state[4*i] = x0[i]
    state[4*i+1] = u0[i]
    state[4*i+2] = y0[i]
    state[4*i+3] = v0[i]

dt = 50
t = np.arange(0, 125000, dt)
# 微分方程组数值解
orbit = integrate.odeint(derivs, state, t)

plt.show()
xMax = np.max(orbit)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(xlim=(-xMax,xMax),ylim=(-xMax,xMax))
ax.grid()

lws = [0.5,2,2,2]
traces = [ax.plot([],[],'-', lw=lws[i])[0] for i in range(len(m))]
pts = [ax.plot([],[], marker='o')[0] for _ in range(len(m))]
time_template = 'time = %.1f d'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

def animate(i):
    for n in range(4):
        traces[n].set_data(orbit[:i,4*n],orbit[:i,4*n+2])
        pts[n].set_data(orbit[i,4*n],orbit[i,4*n+2])
    time_text.set_text(time_template % (i*dt))
    return traces+pts+[time_text]

ani = animation.FuncAnimation(fig, animate, range(len(t)),   
        interval=10, blit=True)
ani.save("tri_5.gif")
plt.show()

接下来还是体验一下行星视角，首先看一下在行星上观察到的恒星们的轨迹

如果看动图可能压迫感会更强一些，这些恒星简直对行星视若无物。

其绘图代码无变化，只需让orbit中的行星位置归零，

for i in range(4):
    orbit[:,4*i] -= orbit[:,0]
    orbit[:,4*i+2] -= orbit[:,2]

由于恒星辐射的功率密度以三次方的形式进行衰减，若假定行星接收到的功率是三颗恒星的叠加，那么就可以画出三体行星所接收的功率变化

由于取了以10为底的对数，所以其峰值功率是最小值的 1 0 7 10^7 107倍，所以这是什么概念呢？

假设在短时间内，功率是均匀的，也就是说单位时间内所爆发出的能量基本是不变的。

一个汉堡的热量大概为2000kJ，那么其 1 0 4 10^4 104倍就是 2 × 1 0 7 k J 2\times10^7kJ 2×107kJ，相当于10发战斧导d。

故而对于三体人而言，严冬之日的一个汉堡包，约等于盛夏之时的十发战斧导d。

所以三体人要是没个超能力什么的，基本上是活不下去的。

P = 0
for i in range(1,4):
    L = np.sqrt(orbit[:,4*i]**2+orbit[:,4*i+2]**2)
    P += 1/L**3

# 为了看上去更清晰，对功率做对数
plt.plot(t,np.log10(P))