基于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


猜你喜欢
- 这篇文章主要介绍了Python爬虫爬取百度搜索内容代码实例,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要
- 想查找mac下python的site-packages目录,如果用网上搜的安装目录之类的去找,进去每一层很多都是软链接形式存在,比较浪费时间
- 这个操作现在看来真没啥难的,但是我找相关的资料真的找了好久。多数大佬都是直接pandas官网甩我脸上,然后举一个入门级的例子。https:/
- 1、权限赋予 说明:mysql部署在服务器A上,内网上主机B通过客户端工具连接服务器A以进行数据库操作,需要服务器A赋予主机B操作mysql
- 为了更好的进行封装,每个实现不同功能的js代码应该有自己的js文件,这样如果一个网页中引用了多个js文件,如下 <script typ
- 本文实例讲述了Python正则表达式实现截取成对括号的方法。分享给大家供大家参考,具体如下:strs = '1(2(3(4(5(67
- 通常我们做统计图的时候需要借助组件来完成例如mschart,aspchart等但是这个类不需要任何组件,而且使用方便,本站测试可用:clsG
- 一般情况下,网站的图片代码是这样的。<img src="./images/test.jpg"
- 一、简介本章内容主要通过具体的简单示例来分析Vue3是如何实现响应式的。理解本章需要了解Vue3的响应式对象。只注重原理设计层面,细节不做太
- 何为样本分布不均:样本分布不均衡就是指样本差异非常大,例如共1000条数据样本的数据集中,其中占有10条样本分类,其特征无论如何你和也无法实
- 1. 视图分离与嵌套在 learnlaravel 文件夹下运行命令:php artisan generate:view admin._lay
- 本文实例讲述了python实现提取百度搜索结果的方法。分享给大家供大家参考。具体实现方法如下:# coding=utf8import url
- 1 create table test(coltest varchar(20))2 实现这一功能 的 sql 语句 s
- By: 吴垠 Date: 2007-09-07 Version: 0.5 Email: lazy.fox.wu#gmail.com Home
- 概述🌱记住日期是有点困难,但我们是程序员,使困难的事情更容易是我们唯一的工作,所以我们不记得日期为什么不自动化这个任务。在这篇文章中,我们将
- 本文实例讲述了C#连接Oracle数据库的方法。分享给大家供大家参考。具体实现方法如下://1、添加引用 System.data.oracl
- 函数描述int(x [,base])将x转换为一个整数long(x [,base] )将x转换为一个长整数float(x)将x转换到一个浮点
- 目录单例模式反射hasattergetattrsetattr总结单例模式一般情况下,类可以生成任意个实例,而单例模式只生成一个实例我们先用单
- 1.file--->settings:2.搜索temlates,点击File and Code Templates,点击右侧的Pyth
- github主页导入:>>> from fuzzywuzzy import fuzz>>> from f