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

三体问题

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

但就拉格朗日方程而言，却并不困难。设 m i , i ∈ { 0 , 1 , 2 , 3 } m_i, i\in\{0,1,2,3\} mi​,i∈{0,1,2,3}表示四颗星体，则对任意星体$i$而言，其动能为

$$T_i=\frac{1}{2}{m}_i({\dot{x}}^2+{\dot{y}}^2)$$

势能为

$$V_i=-\sum _{j\ne i}\frac{G{m}_i{m}_j}{\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}}$$

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

$$\frac{\text{d}}{\text{d}t}\frac{\partial L}{\partial \dot{r}}-\frac{\partial L}{\partial r}=0，r=x_i,y_i$$

有

$$\begin{gathered} {\ddot{x}}_i=-\sum _{i\ne j}\frac{G{m}_j(x_i-x_j)}{{\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}}^3} \\ {\ddot{y}}_i=-\sum _{i\ne j}\frac{G{m}_j(y_i-y_j)}{{\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}}^3} \end{gathered}$$

则令state以如下顺序编排，$s=x_0,{\dot{x}}_0,y_0,{\dot{y}}_0,x_1,{\dot{x}}_1\dots$，则$s_{4i}=x_i,s_{4i+1}={\dot{x}}_i,s_{4i+2}=y_i,s_{4i+3}={\dot{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] = 0dydx[i*4+3] = 0for j in range(N):if i==j:continuedx = state[i*4]-state[j*4]dy = state[i*4+2]-state[j*4+2]L = np.sqrt(dx**2+dy**2)**3dydx[i*4+1] -= G * m[j] * dx / Ldydx[i*4+3] -= G * m[j] * dy / Lreturn dydx

由于三体运动过于放荡不羁，故而随机生成的三体几乎很快就分道扬镳了，所以接下来选择适当位置和重量的三颗恒星。且令万有引力常数以年为时间单位

$$G=4.98\times 10^{-10}k{m}^3{d}^{-2}k{g}^{-1}$$

令恒星质量在$10^{30}kg$的量级，空间距离在$10^{11}km$量级。行星质量在$10^{25}$量级。由于质量相差过大，所以假定行星质量为0也是可以的。

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

首先，星体质量为$m_i$，坐标为$(X_i,Y_i)$，则其重心坐标为

$$x_g=\frac{{\sum }_i{m}_i{X}_i}{\sum m_i}，y_g=\frac{{\sum }_i{m}_i{Y}_i}{\sum m_i}$$

如将坐标系移动到$(x_g,y_g)$，则新坐标系下$x_i=X_i-x_g,y_i=Y_i-x_g$。

则对$m_i$而言，其运动的半径与加速度分别为为

$$\begin{array}{rl}{r}_i& =\sqrt{x_i^2+y_i^2}\\ {\ddot{x}}_i& =\sum _{j\ne i}\frac{G{m}_j(x_j-x_i)}{{\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}}^3}\\ {\ddot{y}}_i& =\sum _{j\ne i}\frac{G{m}_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{array}$$

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

则其初始角度和速度为

$$\begin{array}{rlrl}\cos \theta & =\frac{x_i}{r_i}& & \\ \sin \theta & =\frac{y_i}{r_i}& & \\ u=\dot{x}& =-{\omega }_i{r}_i\sin \theta & =& -{\omega }_i{y}_i\\ v=\dot{y}& ={\omega }_i{r}_i\cos \theta & =& {\omega }_i{x}_i\end{array}$$

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

则其初始化方法为

# 用于初始化星体的质量和位置
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from scipy import integratenp.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 : continuedx = x0[i]-x0[j]dy = y0[i]-y0[j]L = np.sqrt(dx**2+dy**2)**3gm = G * m[j]a[i] += gm * dx / Lb[i] += gm * dy / Lom = 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为底的对数，所以其峰值功率是最小值的$10^7$倍，所以这是什么概念呢？

假设在短时间内，功率是均匀的，也就是说单位时间内所爆发出的能量基本是不变的。一个汉堡的热量大概为2000kJ，那么其$10^4$倍就是$2\times 10^{7}kJ$，相当于10发战斧导弹。

故而对于三体人而言，严冬之日的一个汉堡包，约等于盛夏之时的十发战斧导弹。所以三体人要是没个超能力什么的，基本上是活不下去的。

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))