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

2023-10-29 23:30

本文主要是介绍【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=21mi(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(xixj)2+(yiyj)2 Gmimj

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

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 dtdr˙LrL=0r=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(xixj)2+(yiyj)2 3Gmj(xixj)y¨i=i=j(xixj)2+(yiyj)2 3Gmj(yiyj)

则令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] = 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 × 1 0 − 10 k m 3 d − 2 k g − 1 G=4.98\times10^{-10} km^3d^{-2}kg^{-1} G=4.98×1010km3d2kg1

令恒星质量在 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=miimiXiyg=miimiYi

如将坐标系移动到 ( 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=Xixg,yi=Yixg

则对 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} rix¨iy¨iωi=xi2+yi2 =j=i(xixj)2+(yiyj)2 3Gmj(xjxi)=j=i(xixj)2+(yiyj)2 3Gmj(yjyi)=rix¨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˙=rixi=riyi=ωirisinθ=ωiricosθ==ωiyiωixi

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

在这里插入图片描述

则其初始化方法为

# 用于初始化星体的质量和位置
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为底的对数,所以其峰值功率是最小值的 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发战斧导弹。

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

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

这篇关于【Python】有三颗恒星的三体人很难不产生超能力的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/304085

相关文章

Python爬虫selenium验证之中文识别点选+图片验证码案例(最新推荐)

《Python爬虫selenium验证之中文识别点选+图片验证码案例(最新推荐)》本文介绍了如何使用Python和Selenium结合ddddocr库实现图片验证码的识别和点击功能,感兴趣的朋友一起看... 目录1.获取图片2.目标识别3.背景坐标识别3.1 ddddocr3.2 打码平台4.坐标点击5.图

Python自动化Office文档处理全攻略

《Python自动化Office文档处理全攻略》在日常办公中,处理Word、Excel和PDF等Office文档是再常见不过的任务,手动操作这些文档不仅耗时耗力,还容易出错,幸运的是,Python提供... 目录一、自动化处理Word文档1. 安装python-docx库2. 读取Word文档内容3. 修改

Python重命名文件并移动到对应文件夹

《Python重命名文件并移动到对应文件夹》在日常的文件管理和处理过程中,我们可能会遇到需要将文件整理到不同文件夹中的需求,下面我们就来看看如何使用Python实现重命名文件并移动到对应文件夹吧... 目录检查并删除空文件夹1. 基本需求2. 实现代码解析3. 代码解释4. 代码执行结果5. 总结方法补充在

Python自动化办公之合并多个Excel

《Python自动化办公之合并多个Excel》在日常的办公自动化工作中,尤其是处理大量数据时,合并多个Excel表格是一个常见且繁琐的任务,下面小编就来为大家介绍一下如何使用Python轻松实现合... 目录为什么选择 python 自动化目标使用 Python 合并多个 Excel 文件安装所需库示例代码

Python使用Pandas对比两列数据取最大值的五种方法

《Python使用Pandas对比两列数据取最大值的五种方法》本文主要介绍使用Pandas对比两列数据取最大值的五种方法,包括使用max方法、apply方法结合lambda函数、函数、clip方法、w... 目录引言一、使用max方法二、使用apply方法结合lambda函数三、使用np.maximum函数

Python调用Orator ORM进行数据库操作

《Python调用OratorORM进行数据库操作》OratorORM是一个功能丰富且灵活的PythonORM库,旨在简化数据库操作,它支持多种数据库并提供了简洁且直观的API,下面我们就... 目录Orator ORM 主要特点安装使用示例总结Orator ORM 是一个功能丰富且灵活的 python O

Python使用国内镜像加速pip安装的方法讲解

《Python使用国内镜像加速pip安装的方法讲解》在Python开发中,pip是一个非常重要的工具,用于安装和管理Python的第三方库,然而,在国内使用pip安装依赖时,往往会因为网络问题而导致速... 目录一、pip 工具简介1. 什么是 pip?2. 什么是 -i 参数?二、国内镜像源的选择三、如何

python使用fastapi实现多语言国际化的操作指南

《python使用fastapi实现多语言国际化的操作指南》本文介绍了使用Python和FastAPI实现多语言国际化的操作指南,包括多语言架构技术栈、翻译管理、前端本地化、语言切换机制以及常见陷阱和... 目录多语言国际化实现指南项目多语言架构技术栈目录结构翻译工作流1. 翻译数据存储2. 翻译生成脚本

如何通过Python实现一个消息队列

《如何通过Python实现一个消息队列》这篇文章主要为大家详细介绍了如何通过Python实现一个简单的消息队列,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录如何通过 python 实现消息队列如何把 http 请求放在队列中执行1. 使用 queue.Queue 和 reque

Python如何实现PDF隐私信息检测

《Python如何实现PDF隐私信息检测》随着越来越多的个人信息以电子形式存储和传输,确保这些信息的安全至关重要,本文将介绍如何使用Python检测PDF文件中的隐私信息,需要的可以参考下... 目录项目背景技术栈代码解析功能说明运行结php果在当今,数据隐私保护变得尤为重要。随着越来越多的个人信息以电子形