从零开始数值模拟三体问题||Python(For新手)

//最近试着自学了python对初值问题的数值模拟。

//原理也非常简单，把微分方程近似为差分方程，然后用循环语句计算质点位置的变化。

//本程序基于python的xlrd，numpy，matplotlib，mpl_toolkits.mplot3D库实现。

//xlrd可以从excel导入数据，numpy处理向量运算，matplotlib和mpl_toolkits绘制图像。

0. 理论基础

三体系统是著名的混沌系统，三个质点仅受万有引力作用而产生复杂的运动，且决定系统运动的微分方程无法求得解析解，只能数值求解。同时，系统的后续运动对初值高度敏感，随着时间推移数值近似的解误差也越来越大。虽然对混沌系统的详细研究相当复杂，有包括庞加莱截面法、功率谱分析法、分形维数、李雅普诺夫指数等刻画混沌现象的方法，但不影响我们基于基本的微分方程对该系统进行简单的数值模拟。

设三个星体质量分别为，位置用矢量表示。三体系统的微分方程如下：

在现实中万有引力属于相对微弱的作用力(相对其他基本力)，因此引力的作用也通常是质量极大的天体在相当长的时间里相互作用得以体现(相对地球上的日常而言)。所以这个模型中我们简单起见直接令，所有量纲相关的讨论都不考虑。

作为数值模拟的需要，我们将微分方程化为差分方程，以1质点为例：

受力

变速

位移

其中是程序中设定的一个很小的时间间隔。越小则误差越小，但计算量越大。

I. 准备工作

开始之前你需要确认你的python装有前面提到的xlrd，numpy，matplotlib，mpl_toolkits几个库。如果你使用的IDE是Anaconda，这些包都是自带的，不需要额外安装。关于其他如Pycharm, Visual Studio等环境我也不是很了解，但总之如果你缺少某个包(例如numpy)，只要在命令行中运行以下命令：

系统就会自动从网络上搜索相关库并安装。确认所有需要的库都安装后就可以开始了！

II. 开始编程

①包与数据导入

我们首先在程序中导入所有需要的库：

由于本系统需要的初始条件高达21个(每个质点的质量、初始位置三个分量、初速度三个分量)，我们如果利用input()命令在程序运行时逐个输入麻烦且不直观，我们选择把初值写在excel文件中并利用xlrd包导入。首先找一个地方创建一个Excel文件DataInput.xlsx，例如在D盘的文档中。Python文件中利用以下命令：

就可以将数据来源指定为DataInput.xlsx的sheet1。在excel中，按如图方式将三个质点的质量、初位置、初速度写入，保存。(你可以根据自己的喜好指定初始值，图中的数据仅供参考)

接下来我们将这些数据存储在变量m, 数组r,v中，作为后续计算的准备工作：

注意虽然我们的数据是从第2行第2列开始的，但代码中是从(1,1)开始导入数据的，这是因为大部分编程语言中数组或类似结构下标都默认从0开始，下标要注意减1.

② 计算准备

我们需要准备九个空列表存放各个质点坐标随时间的变化：

设置和总步数：

创建一个三维图，用于可视化结果：

我们定义几个函数，求点间的施力：

四个函数的功能分别是：输入两个位置返回距离；输入距离和质量返回引力大小；输入两个位置返回连接它们的单位矢量(以数组形式)；调用引力大小和单位矢量函数返回引力矢量(数组形式)。其中，调用了numpy中的数学、数组运算函数。

接下来，定义两个void函数，用于进行质点的移动和受力改变速度：

③ 计算

准备工作都做好了，现在只需要一次循环语句完成所有计算即可。

④ 绘图

所有数据的计算都完成了，接下来就简单了，以下命令将几个质点的轨迹绘制在前面定义的三维图中：

最后输入预览命令

我们就会成功看到前面的封面图片了！

可以看到，确实很混沌。

如果你希望保存这个图片，即可以直接在预览窗口上点击保存，也可以用命令

来保存这个结果。

那么以上就是这样一个新手向的Python数值计算笔记了。当然本人也只是Python初学者，偶尔用它来数据分析，这个程序可能没有写出最简洁高效的代码，如果你有更好的方法也欢迎指出。