脚本之家 服务器常用软件
快捷导航
您的位置：首页 → 脚本专栏 → python → Python三体运动

基于Python实现模拟三体运动的示例代码

 更新时间：2023年03月10日 08:44:01   作者：微小冷  
此前所做的一切三体和太阳系的动画，都是基于牛顿力学的，而且直接对微分进行差分化，从而精度非常感人，用不了几年就得撞一起去。所以本文来用Python重新模拟一下三体运动，感兴趣的可以了解一下

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

拉格朗日方程

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

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

首先假设三个质点的质量分别为m1, m2,m3,坐标为x→1​,x→2​,x→3,​质点速度可以表示为x → ˙.假设三体在二维平面上运动，则第i个质点的动能为

引力势能为

其中G为万有引力常量，rij为质点i,j之间的距离，则系统的拉格朗日量为

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

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

推导方程组

对于三体系统而言，总计有3个粒子，每个粒子有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，其实质是系统的动能减去势能，对于上面构建的三体系统而言，动能和势能可分别表示为

计算每个质点的动能和势能。动能是由速度决定的，而速度是由位置对时间的导数决定的。我们可以用 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]

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

三个粒子的每一个坐标维度，都可以列出一组拉格朗日方程，所以总共有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]

记xij=xi−xj,yij=yi−yj ，则

微分方程算法化

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

为此，需要将上述方程组再行拆分，以消去其中的二次导数，以x1为例，令u1=dx1/dt ，则此方程变为方程组

由于三体系统中有3个粒子，共6个独立变量，所以要列12个方程。记

则odeint输入的y的形式为

从而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')

到此这篇关于基于Python实现模拟三体运动的示例代码的文章就介绍到这了,更多相关Python三体运动内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！

您可能感兴趣的文章:

相关文章

  • Python异常对象Exception基础类异常捕捉

    Python异常对象Exception基础类异常捕捉

    这篇文章主要为大家介绍了Python异常对象异常捕捉及Exception基础类，有需要的朋友可以借鉴参考下，希望能够有所帮助，祝大家多多进步，早日升职加薪
    2022-06-06
  • python爬虫之你好,李焕英电影票房数据分析

    python爬虫之你好,李焕英电影票房数据分析

    这篇文章主要介绍了python爬虫之你好,李焕英电影票房数据分析,文中有非常详细的代码示例,对正在学习python爬虫的小伙伴们有一定的帮助,需要的朋友可以参考下
    2021-04-04
  • Python批量发送post请求的实现代码

    Python批量发送post请求的实现代码

    昨天学了一天的Python（我的生产语言是java，也可以写一些shell脚本，算有一点点基础），今天有一个应用场景，就正好练手了
    2018-05-05
  • 关于Python中几种队列Queue用法区别

    关于Python中几种队列Queue用法区别

    这篇文章主要介绍了关于Python中几种队列Queue用法区别,queue队列中的put()或者get()方法中都提供了timeout参数，利用这个参数可以有效解决上述消除不能消费和线程一直阻塞问题,需要的朋友可以参考下
    2023-05-05
  • Python实现全自动输入文本的示例详解

    Python实现全自动输入文本的示例详解

    这篇文章主要和大家分享一个Python全自动输入文本的脚本，可以实现自动用Notepad++打开文本文件，然后自动输入文本，最后保存并关闭文件，从而实现全面自动化处理文本，希望对大家有所帮助
    2022-11-11
  • Python实现在PDF中插入单图像水印和平铺图像水印

    Python实现在PDF中插入单图像水印和平铺图像水印

    这篇文章主要为大家详细介绍了如何使用Python实现在PDF中插入单图像水印和平铺图像水印,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下
    2024-04-04
  • Python 合并多个TXT文件并统计词频的实现

    Python 合并多个TXT文件并统计词频的实现

    这篇文章主要介绍了Python 合并多个TXT文件并统计词频的实现，文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值，需要的朋友们下面随着小编来一起学习学习吧
    2019-08-08
  • python数据结构的排序算法

    python数据结构的排序算法

    下面是是对python数据结构的排序算法的一些讲解及示意图，感兴趣的小伙伴一起来学习吧
    2021-08-08
  • PyQt5之如何设置QWidget窗口背景图片问题

    PyQt5之如何设置QWidget窗口背景图片问题

    这篇文章主要介绍了PyQt5之如何设置QWidget窗口背景图片问题，具有很好的参考价值，希望对大家有所帮助。如有错误或未考虑完全的地方，望不吝赐教
    2023-06-06
  • Python ChineseCalendar包主要类和方法详解

    Python ChineseCalendar包主要类和方法详解

    ChineseCalendar 是一个 Python 包，用于获取中国传统日历信息。这个包提供了中国农历、二十四节气、传统节日、黄历等信息，这篇文章主要介绍了Python ChineseCalendar包简介,需要的朋友可以参考下
    2023-03-03

最新评论