基于Python实现模拟三体运动的示例代码
作者:微小冷 发布时间:2022-03-29 21:40:37
温馨提示,只想看图的画直接跳到最后一节
拉格朗日方程
此前所做的一切三体和太阳系的动画,都是基于牛顿力学的,而且直接对微分进行差分化,从而精度非常感人,用不了几年就得撞一起去。
为了给三体人提供一个更加有价值的推导,这次通过求解拉格朗日方程的数值解来实现。
首先假设三个质点的质量分别为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')
来源:https://blog.csdn.net/m0_37816922/article/details/129430989
猜你喜欢
- 一、前言我们一般在做接口关联时,会通过保存中间变量实现接口关联,在关联时就需要用到变量提取,那今天我们就介绍接口自动化中变量提取的两大神器:
- 我就废话不多说了,大家还是直接看代码吧~func ReadLine(fileName string) ([]string,error){f,
- 本文为大家分享了Python2.7与Python3.6环境切换的具体方法,供大家参考,具体内容如下系统支持为:Ubuntu18.04系统默认
- 很多小伙伴都会有这样的问题,说一个ip地址十分钟内之内注册一次,用来防止用户来重复注册带来不必要的麻烦逻辑:取ip,在数据库找ip是否存在,
- 很多书籍里面讲的Python备份都是在linux下的,而在xp上测试一下也可以执行备份功能,代码都差不多相同,就是到执行打包的时候是不一样的
- 大家知道,在js里encodeURIComponent 方法是一个比较常用的编码方法,但因工作需要,在asp里需用到此方法,查了好多资料,没
- 引入在通过scrapy框架进行某些网站数据爬取的时候,往往会碰到页面动态数据加载的情况发生,如果直接使用scrapy对其url发请求,是绝对
- 印刷和网络是不一样的。传统的布局排版并不适于网络,因为传统的印刷布局,几乎只想要什么样的平面效果都能很好的达到,但在网络上设计就很困难,尽管
- 在采集美女站时,需要对关键词进行分词,最终采用的是python的结巴分词方法。中文分词是中文文本处理的一个基础性工作,结巴分词利用进行中文分
- 本文实例讲述了django实现分页的方法。分享给大家供大家参考。具体如下:Python代码如下:#!/usr/bin/env python#
- 1. __init__ 初始化文件路径,关键字1,关键字2;2. key_match 使用with open 方法,以二进制方式(也可以改成
- 利用Python写了一个小脚本想要传给使用Windows但没有装Python的朋友执行,这时候就可以利用将档案包装成exe档案,让没有Pyt
- 用python进行线性回归分析非常方便,有现成的库可以使用比如:numpy.linalog.lstsq例子、scipy.stat
- requests接口测试的介绍requests是一个很实用的Python HTTP客户端库,编写爬虫和测试服务器响应数据时经常会用到,Req
- 前言Tkinter(即 tk interface) 是 Python 标准 GUI 库,简称 “Tk&rdquo
- scrapy框架只能爬取静态网站。如需爬取 * 站,需要结合着selenium进行js的渲染,才能获取到动态加载的数据。如何通过seleni
- 举例: 340%60 = 40 ,怎么来的?340 - 60*5 = 40340 - (比340小的那个可以被60整除的正整数) =. 40
- 问题:关于如何生成随机记录(二)如何从指定表中随机抽取一定量的记录?sql server 中 select top 10 * fr
- 本文实例为大家分享了python修改装饰器中参数的具体代码,供大家参考,具体内容如下案例: &
- 在这之前我们先回顾以前用php导出excel,我直接写成方法在这里:public static function phpExcelList(