三体到底是啥?用Python跑一遍就明白了

文章目录

    • 拉格朗日方程
    • 推导方程组
    • 微分方程算法化
    • 求解+画图
    • 动图绘制

温馨提示,只想看图的画直接跳到最后一节

拉格朗日方程

此前所做的一切三体和太阳系的动画,都是基于牛顿力学的,而且直接对微分进行差分化,从而精度非常感人,用不了几年就得撞一起去。

为了给三体人提供一个更加有价值的推导,这次通过求解拉格朗日方程的数值解来实现。

首先假设三个质点的质量分别为 m 1 , m 2 , m 3 m_1, m_2, m_3 m1,m2,m3,坐标为 x ⃗ 1 , x ⃗ 2 , x ⃗ 3 \vec x_1, \vec x_2, \vec x_3 x 1,x 2,x 3,质点速度可以表示为 x ⃗ ˙ \dot{\vec x} x ˙。假设三体在二维平面上运动,则第 i i i个质点的动能为

T i = 1 2 m i ( x ˙ i 2 + y ˙ i 2 ) T_i=\frac{1}{2}m_i(\dot x_i^2+\dot y_i^2) Ti=21mi(x˙i2+y˙i2)

引力势能为 − G m i m j r i j -G\frac{m_im_j}{r_{ij}} Grijmimj,其中 G G G为万有引力常量, r i j r_{ij} rij为质点 i , j i,j i,j之间的距离,则系统的拉格朗日量为

L = ∑ i 1 2 m i ( x ˙ i 2 + y ˙ i 2 ) − ∑ i ≠ j G m i m j ∥ x ⃗ i − x ⃗ j ∥ L=\sum_i\frac{1}{2}m_i(\dot x_i^2+\dot y_i^2)-\sum_{i\not=j}G\frac{m_im_j}{\Vert\vec x_i-\vec x_j\Vert} L=i21mi(x˙i2+y˙i2)i=jGx ix jmimj

有了拉格朗日量,将其带入拉格朗日方程

d d t ∂ L ∂ x ˙ i − ∂ L ∂ x i = 0 \frac{\text d}{\text dt}\frac{\partial L}{\partial\dot x_i}-\frac{\partial L}{\partial x_i}=0 dtdx˙iLxiL=0

就可以得到拉格朗日方程组。

推导方程组

对于三体系统而言,总计有3个粒子,每个粒子有 x , y x,y x,y两个自由度,也就是说最后会得到6组方程。考虑到公式推导过程中可能会出现错误,所以下面采用sympy来进行公式推导。

首先定义符号变量

from sympy import symbols
from sympy.physics.mechanics import dynamicsymbols
m = symbols('m1:4')
x = dynamicsymbols('x1:4')
y = dynamicsymbols('y1:4')

接下来,需要构造系统的拉格朗日量 L L L,其实质是系统的动能减去势能,对于上面构建的三体系统而言,动能和势能可分别表示为

计算每个质点的动能和势能。动能是由速度决定的,而速度是由位置对时间的导数决定的。我们可以用 sympy 的 diff 函数来求导:

from sympy import diff
# 此为速度的平方
v2 = [diff(x[i],t)**2 + diff(y[i])**2 for i in range(3)]
T = 0
for i in range(3):
    T += m[i]*v2[i]/2

势能是由万有引力决定的,而万有引力是由两个质点之间的距离决定的。我们可以用 sympy 的 sqrt 函数来求距离:

from sympy import sqrt,cos
G = symbols('G') # 引力常数
ijs = [(0,1), (0,2),(1,2)]
dij = [sqrt((x[i]-x[j])**2+(y[i]-y[j])**2) for i,j in ijs]
U = 0
for k in range(3):
    i,j = ijs[k]
    U -= G*m[i]*m[j]/dij[k]

有了动能和势能,就可以愉快地求拉格朗日量了,有了拉格朗日量,就可以列拉格朗日方程了

L = T − U d L d x i − d d t ∂ L ∂ x ˙ i L = T - U\\ \frac{\text dL}{\text dx_i}-\frac{\text d}{\text dt}\frac{\partial L}{\partial \dot x_i} L=TUdxidLdtdx˙iL

三个粒子的每一个坐标维度,都可以列出一组拉格朗日方程,所以总共有6个拉格朗日方程组

from sympy import solve
L = T - U
eqLag = lambda x : diff(L, x)-diff(diff(L, diff(x, t)), t)
# 拉格朗日方程组
eqs = [eqLag(xi) for xi in x+y]

x i j = x i − x j , y i j = y i − y j x_{ij}=x_i-x_j, y_{ij}=y_i-y_j xij=xixj,yij=yiyj,则

− G m 1 m 2 x 12 ( x 12 2 + y 12 2 ) 3 2 + − G m 1 m 3 x 13 ( x 13 2 + y 13 2 ) 3 2 − m 1 d 2 d t 2 x 1 = 0 G m 1 m 2 x 12 ( x 12 2 + y 12 2 ) 3 2 + − G m 2 m 3 x 23 ( x 23 2 + y 23 2 ) 3 2 − m 2 d 2 d t 2 x 2 = 0 G m 1 m 3 x 13 ( x 13 2 + y 13 2 ) 3 2 + G m 2 m 3 x 23 ( x 23 2 + y 23 2 ) 3 2 − m 3 d 2 d t 2 x 3 = 0 − G m 1 m 2 y 12 ( x 12 2 + y 12 2 ) 3 2 + − G m 1 m 3 y 13 ( x 13 2 + y 13 2 ) 3 2 − m 1 d 2 d t 2 y 1 = 0 G m 1 m 2 y 12 ( x 12 2 + y 12 2 ) 3 2 + − G m 2 m 3 y 23 ( x 23 2 + y 23 2 ) 3 2 − m 2 d 2 d t 2 y 2 = 0 G m 1 m 3 y 13 ( x 13 2 + y 13 2 ) 3 2 + G m 2 m 3 y 23 ( x 23 2 + y 23 2 ) 3 2 − m 3 d 2 d t 2 y 3 = 0 \frac{-G m_1 m_2x_{12}}{\left(x_{12}^{2} + y_{12}^{2}\right)^{\frac{3}{2}}} + \frac{-G m_1 m_{3}x_{13}}{\left(x_{13}^{2} + y_{13}^{2}\right)^{\frac{3}{2}}} - m_1 \frac{d^{2}}{d t^2} x_1=0\\ \frac{G m_1 m_2 x_{12}}{\left(x_{12}^{2} + y_{12}^{2}\right)^{\frac{3}{2}}} + \frac{-G m_2 m_{3}x_{23}}{\left(x_{23}^{2} + y_{23}^{2}\right)^{\frac{3}{2}}} - m_2 \frac{d^{2}}{d t^2} x_2=0\\ \frac{G m_1 m_{3} x_{13}}{\left(x_{13}^{2} + y_{13}^{2}\right)^{\frac{3}{2}}} + \frac{G m_2 m_{3} x_{23}}{\left(x_{23}^{2} + y_{23}^{2}\right)^{\frac{3}{2}}} - m_{3} \frac{d^{2}}{d t^2} x_{3}=0\\ \frac{-G m_1 m_2 y_{12}}{\left(x_{12}^{2} + y_{12}^{2}\right)^{\frac{3}{2}}} + \frac{-G m_1 m_{3} y_{13}}{\left(x_{13}^{2} + y_{13}^{2}\right)^{\frac{3}{2}}} - m_1 \frac{d^{2}}{d t^2} y_1=0\\ \frac{G m_1 m_2 y_{12}}{\left(x_{12}^{2} + y_{12}^{2}\right)^{\frac{3}{2}}} + \frac{-G m_2 m_{3}y_{23}}{\left(x_{23}^{2} + y_{23}^{2}\right)^{\frac{3}{2}}} - m_2 \frac{d^{2}}{d t^2} y_2=0\\ \frac{G m_1 m_{3} y_{13}}{\left(x_{13}^{2} + y_{13}^{2}\right)^{\frac{3}{2}}} + \frac{G m_2 m_{3} y_{23}}{\left(x_{23}^{2} + y_{23}^{2}\right)^{\frac{3}{2}}} - m_{3} \frac{d^{2}}{d t^2} y_{3}=0\\ (x122+y122)23Gm1m2x12+(x132+y132)23Gm1m3x13m1dt2d2x1=0(x122+y122)23Gm1m2x12+(x232+y232)23Gm2m3x23m2dt2d2x2=0(x132+y132)23Gm1m3x13+(x232+y232)23Gm2m3x23m3dt2d2x3=0(x122+y122)23Gm1m2y12+(x132+y132)23Gm1m3y13m1dt2d2y1=0(x122+y122)23Gm1m2y12+(x232+y232)23Gm2m3y23m2dt2d2y2=0(x132+y132)23Gm1m3y13+(x232+y232)23Gm2m3y23m3dt2d2y3=0

微分方程算法化

接下来就要调用Python的odeint来计算这个微分方程组的数值解,odeint的调用方法大致为odeint(func, y, t, args),其中func是一个函数,这个函数必须为func(y,t,...),且返回值为 d y d t \frac{\text dy}{\text dt} dtdy

为此,需要将上述方程组再行拆分,以消去其中的二次导数,以 x 1 x_1 x1为例,令 u 1 = d x 1 d t u_1=\frac{\text dx_1}{\text dt} u1=dtdx1,则此方程变为方程组

x ˙ 1 ( t ) = u 1 ( t ) u ˙ 1 ( t ) = − G m 1 m 2 x 12 ( x 12 2 + y 12 2 ) 3 2 + − G m 1 m 3 x 13 ( x 13 2 + y 13 2 ) 3 2 \begin{aligned} \dot x_1(t)&=u_1(t)\\ \dot u_1(t)&= \frac{-G m_1 m_2x_{12}}{\left(x_{12}^{2} + y_{12}^{2}\right)^{\frac{3}{2}}} + \frac{-G m_1 m_{3}x_{13}}{\left(x_{13}^{2} + y_{13}^{2}\right)^{\frac{3}{2}}}\\ \end{aligned} x˙1(t)u˙1(t)=u1(t)=(x122+y122)23Gm1m2x12+(x132+y132)23Gm1m3x13

由于三体系统中有3个粒子,共6个独立变量,所以要列12个方程。记 u ( t ) = t e x t d x d t , v ( t ) = d y d t u(t)=\frac{text dx}{\text dt}, v(t)=\frac{\text dy}{\text dt} u(t)=dttextdx,v(t)=dtdy,则odeint输入的y的形式为

[ x 1 , x 2 , x 3 , y 1 , y 2 , y 3 , u 1 , u 2 , u 3 , v 1 , v 2 , v 3 ] [x_1, x_2, x_3, y_1, y_2, y_3, u_1, u_2, u_3, v_1, v_2, v_3] [x1,x2,x3,y1,y2,y3,u1,u2,u3,v1,v2,v3]

从而func的具体形式为

import numpy as np
dxy = lambda x,y : np.sqrt(x**2+y**2)**(3/2)
def triSys(Y, t, m, G):
    jk = [(1,2),(0,2),(0,1)]
    x,y = Y[:3], Y[3:6]
    u,v = Y[6:9], Y[9:]
    du, dv = [], []
    for i in range(3):
        j, k = jk[i]
        xji, xki = x[j]-x[i], x[k]-x[i]
        yji, yki = y[j]-y[i], y[k]-y[i]
        dji, dki = dxy(xji, yji), dxy(yji, yki)
        mji, mki = G*m[i]*m[j], G*m[i]*m[k]
        du.append(mji*xji/dji + mki*xki/dki)
        dv.append(mji*yji/dji + mki*yki/dki)
    dydt = [*u, *v, *du, *dv]
    return dydt

求解+画图

接下来就是见证奇迹的时刻,首先创建一个随机的起点,作为三体运动的初值,然后带入开整就完事儿了

from scipy.integrate import odeint
np.random.seed(42)
y0 = np.random.rand(12)
m = np.random.rand(3)
t = np.linspace(0, 20, 1001)
sol = odeint(triSys, y0, t, args=(m, 1))

然后绘制一下这三颗星的轨迹

import matplotlib.pyplot as plt
plt.plot(sol[:,0], sol[:,3])
plt.plot(sol[:,1], sol[:,4])
plt.plot(sol[:,2], sol[:,5])
plt.show()

在这里插入图片描述

光是看这个轨迹就十分惊险了有木有。

如果把其中的第一颗星作为坐标原点,那么另外两颗星的轨迹大致为

plt.plot(sol[:,1]-sol[:,0], sol[:,4]-sol[:,3])
plt.plot(sol[:,2]-sol[:,0], sol[:,5]-sol[:,3])
plt.scatter([0],[0], c='g', marker='*')
plt.show()

结果为

在这里插入图片描述

动图绘制

最后,以中间这颗星为原点,绘制一下另外两颗星运动的动态过程

import matplotlib.animation as animation 

fig = plt.figure(figsize=(9,4))
ax = fig.add_subplot(xlim=(-1.8,1.8),ylim=(-1.8,1.5))
ax.grid()

traces = [ax.plot([],[],'-',lw=0.5)[0] for _ in range(2)]
pts = [ax.plot([],[] ,marker='*')[0] for _ in range(2)]
ax.plot([0],[0], marker="*", c='r')

X1 = sol[:,1]-sol[:,0]
Y1 = sol[:,4]-sol[:,3]
X2 = sol[:,2]-sol[:,0]
Y2 = sol[:,5]-sol[:,3]

def animate(n):
    traces[0].set_data(X1[:n], Y1[:n])
    traces[1].set_data(X2[:n], Y2[:n])
    pts[0].set_data([X1[n], Y1[n]])
    pts[1].set_data([X2[n], Y2[n]])
    return traces + pts

ani = animation.FuncAnimation(fig, animate, 
    range(1000), interval=10, blit=True)
ani.save('tri.gif')

在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mfbz.cn/a/335.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

如何用Python求解微分方程组

文章目录odeint简介示例odeint简介 scipy文档中将odeint函数和ode, comples_ode这两个类称为旧API,是scipy早期使用的微分方程求解器,但由于是Fortran实现的,尽管使用起来并不方便,但速度没得说,所以有的时候还挺推荐…

Vite4 + Vue3 + vue-router4 动态路由

动态路由,基本上每一个项目都能接触到这个东西,通俗一点就是我们的菜单是根据后端接口返回的数据进行动态生成的。表面上是对菜单的一个展现处理,其实内部就是对router的一个数据处理。这样就可以根据角色权限或者一些业务上的需求&#xff0…

机器学习入门——线性回归

线性回归什么是线性回归?回归分析:线性回归:回归问题求解单因子线性回归简单实例评估模型表现可视化模型展示多因子线性回归什么是线性回归? 回归分析: 根据数据,确定两种或两种以上变量间相互依赖的定量…

自学大数据第六天~HDFS命令(一)

HDFS常用命令 查看hadoop版本 version hadoop version注意,没有 ‘-’ [hadoopmaster ~]$ hadoop version Hadoop 3.3.4 Source code repository https://github.com/apache/hadoop.git -r a585a73c3e02ac62350c136643a5e7f6095a3dbb Compiled by stevel on 2022-07-29T12:3…

【电赛MSP430系列】GPIO、LED、按键、时钟、中断、串口、定时器、PWM、ADC

文章目录MSP430一、GPIO二、点亮LED三、按键控制LED四、更改主时钟五、串口通信六、串口中断七、外部中断八、定时器九、定时器中断十、PWM十一、ADCMSP430 MSP430 是德州仪器(TI)一款性能卓越的超低功耗 16 位单片机,自问世以来&#xff0c…

程序员的逆向思维

前要: 为什么你读不懂面试官提问的真实意图,导致很难把问题回答到面试官心坎上? 为什么在面试结束时,你只知道问薪资待遇,不知道如何高质量反问? 作为一名程序员,思维和技能是我们职场生涯中最重要的两个方面。有时候…

【微信小程序】-- 网络数据请求(十九)

💌 所属专栏:【微信小程序开发教程】 😀 作  者:我是夜阑的狗🐶 🚀 个人简介:一个正在努力学技术的CV工程师,专注基础和实战分享 ,欢迎咨询! &…

到底什么是跨域,如何解决跨域(常见的几种跨域解决方案)?

文章目录1、什么是跨域2、解决跨域的几种方案2.1、JSONP 方式解决跨域2.2、CORS 方式解决跨域(常见,通常仅需服务端修改即可)2.3、Nginx 反向代理解决跨域(推荐使用,配置简单)2.4、WebSocket 解决跨域2.5、…

软测面试了一个00后,绝对能称为是内卷届的天花板

前言 公司前段缺人,也面了不少测试,结果竟然没有一个合适的。一开始瞄准的就是中级的水准,也没指望来大牛,提供的薪资也不低,面试的人很多,但平均水平很让人失望。令我印象最深的是一个00后测试员&#xf…

【JavaScript 逆向】百度旋转验证码逆向分析

声明本文章中所有内容仅供学习交流,相关链接做了脱敏处理,若有侵权,请联系我立即删除!案例目标爱企查百度安全验证百度搜索:aHR0cHM6Ly93YXBwYXNzLmJhaWR1LmNvbS9zdGF0aWMvY2FwdGNoYS8以上均做了脱敏处理,B…

操作系统(2.2)--进程的描述与控制

目录 二、进程的描述 1.进程的定义和特征 1.1进程的定义 1.2进程的特征 2.进程的基本状态及转换 2.1进程的三种基本状态 2.2 三种基本状态的转换 2.3创建状态和中止状态 3.挂起操作和进程状态的转换 3.1 挂起状态的引入 3.2 引入挂起操作后三个进程状态的转换 …

07从零开始学Java之如何正确的编写Java代码?

作者:孙玉昌,昵称【一一哥】,另外【壹壹哥】也是我哦CSDN博客专家、万粉博主、阿里云专家博主、掘金优质作者前言在上一篇文章中,壹哥带领大家开始编写了第一个Java案例,在我们的cmd命令窗口中输出了”Hello World“这…

【蓝桥杯-筑基篇】常用API 运用(1)

🍓系列专栏:蓝桥杯 🍉个人主页:个人主页 目录 🍍1.输入身份证,判断性别🍍 🍍2.输入英语句子,统计单词个数🍍 🥝3.加密解密🥝 🌎4.相邻重复子串…

【6G 新技术】6G数据面介绍

博主未授权任何人或组织机构转载博主任何原创文章,感谢各位对原创的支持! 博主链接 本人就职于国际知名终端厂商,负责modem芯片研发。 在5G早期负责终端数据业务层、核心网相关的开发工作,目前牵头6G算力网络技术标准研究。 博客…

订单30分钟未支付自动取消怎么实现?

目录了解需求方案 1:数据库轮询方案 2:JDK 的延迟队列方案 3:时间轮算法方案 4:redis 缓存方案 5:使用消息队列了解需求在开发中,往往会遇到一些关于延时任务的需求。例如生成订单 30 分钟未支付&#xff0…

HashMap的实际开发使用

目 录 前言 一、HashMap是什么? 二、使用步骤 1.解析一下它实现的原理 ​编辑 2.实际开发使用 总结 前言 本章,只是大概记录一下hashMap的简单使用方法,以及理清一下hashMap的put方法的原理,以及get方法的原理。 一、Has…

如何使用 Python 检测和识别车牌(附 Python 代码)

文章目录创建Python环境如何在您的计算机上安装Tesseract OCR?技术提升磨砺您的Python技能车牌检测与识别技术用途广泛,可以用于道路系统、无票停车场、车辆门禁等。这项技术结合了计算机视觉和人工智能。 本文将使用Python创建一个车牌检测和识别程序。…

【Linux】目录和文件的权限

Linux中的权限有什么作用Linux权限管理文件访问者的分类文件类型和访问权限(事物属性)**文件权限值的表示方法**文件访问权限的相关设置方法chmodchownchgrpumaskumask使用 sudo分配权限目录的权限Linux中的权限有什么作用 Linux下有两种用户&#xff1…

【C缺陷与陷阱】----语义“陷阱”

💯💯💯 本篇处理的是有关语义误解的问题:即程序员的本意是希望表示某种事物,而实际表示的却是另外一种事物。在本篇我们假定程序员对词法细节和语法细节的理解没有问题,因此着重讨论语义细节。导言&#xf…

SignalR+WebRTC技术实现音视频即时通讯功能

一、建立信令服务器 1、后台项目中新建一个对应的集线器类,取名VideoHub,并继承Hub类,Hub是SignalR的一个组件,它使用RPC接收从客户端发送来的消息,也能把消息发送给客户端。 2、VideoHub中定义一个静态的Concurrent…
最新文章