The N-Body Problem
三体看了这么多年了，自然也会对三体问题产生兴趣。虽然大刘在书里面说三体问题无解（我被骗了很多年），但本着科学的实证精神（在知乎上刷到了一篇关于三体问题的科普），我去系统的调查了一番，发现三体问题实际上是有解析解的，而且庞加莱老早以前就证明了，只不过解析解是级数形式的，而且通常收敛性很差，基本没法用，所以目前的主流方法还是数值模拟。
两个月前，体校有大佬打算用Java写黑森的多人联网模拟程序（虽然现在貌似遇到了点小问题），中间顺带写了个多体模拟程序，可以做到三千体量级的模拟，我也因此起了这方面的想法。
我正式开始写程序是在七月初，用的是Mathematica，暂时命名为N-Body，到目前为止经过了多次改进，差不多到我的瓶颈了，因此把这部分的经历记录一下，一方面作为科普，另一方面也希望有大佬能够提出改进的意见，使程序继续完善
首先要从基本原理说起，我们知道，对于宏观低速的天体系统，最适合的描述方法是万有引力定律与牛顿运动定律。假定一个封闭的多体引力系统，其中的所有天体只受其他天体的引力作用，由引力决定天体此时的加速度，由加速度决定天体的速度演化，由速度决定天体的运动轨迹。
用数学语言来表述就是（以三体系统为例）：
因为里面的r1、r2、r3都是三维向量，因此这个方程组可以被分解到三个方向上，在每个方向都满足以上的形式。
因为里面的r1、r2、r3都是三维向量，因此这个方程组可以被分解到三个方向上，在每个方向都满足以上的形式。
目前网上的大部分三体模拟代码都是基于这一点，用ODE库或者内置的函数直接求该方程组的解数值，然后取点画图。这种方法的优点是代码简单而且精度高，缺点是不能随意地扩展体系大小，如果你要模拟3000体的运动，那么就需要解组包含3000个方程的方程组，而且每一个方程都有2999项
所以，为了实现体系的任意扩展，我们必须脱离方程组的形式，从纯数值的角度来解决这个问题。
因为我此前也没有相关经验，所以为了方便上手，我最初采用的方法是欧拉法，也就是：
因为我此前也没有相关经验，所以为了方便上手，我最初采用的方法是欧拉法，也就是：
程序里的数据是以数组的形式存在的（mathematica里的数组由{}分隔），其中质量是一维数组，它的长度决定了体系所包含的天体数目，坐标和速度则是三维数组，三维数组第一层的元素表示某一时刻所有天体的坐标或速度。
程序能够实现体系的扩展主要依靠加速度计算函数f，这个函数的输入和输出都是二维数组（也就是矩阵），在这里，输入的二维数组代表各个天体的坐标，输出的二维数组则代表各个天体的加速度
mathematica所使用的wolfram语言是一种逐行解释型语言，语法非常贴近自然语言，但它的底层语言是c。通过调用超过6000个build-in函数，mathematica可以用非常简单的代码实现许多复杂的功能，而且拥有不错的速度。
我在构造f函数时就使用了六个内置函数：Compile,Lenth,Table,Total,If,Norm
我在构造f函数时就使用了六个内置函数：Compile,Lenth,Table,Total,If,Norm
Compile函数的作用是预先编译自定义函数，得到一个运行在底层逻辑上模块，节约函数的运算时间。
Lenth和If的作用很简单，一个是计算数组第一层的长度，另一个就是基础的逻辑判断，条件为真执行条件后第一个语句，否则执行第二个。
Table是mma里面使用十分广泛的一个函数，虽然翻译成表格，但更合适的解释是“迭代器”，它在mma里的地位相当于c、py和Java里的for循环。
Total可以理解为∑，实际上，mma在处理符号运算时就是把∑作为Total处理。
Norm的作用则是得到输入向量的模长。
Lenth和If的作用很简单，一个是计算数组第一层的长度，另一个就是基础的逻辑判断，条件为真执行条件后第一个语句，否则执行第二个。
Table是mma里面使用十分广泛的一个函数，虽然翻译成表格，但更合适的解释是“迭代器”，它在mma里的地位相当于c、py和Java里的for循环。
Total可以理解为∑，实际上，mma在处理符号运算时就是把∑作为Total处理。
Norm的作用则是得到输入向量的模长。
综上所述，f函数的逻辑过程是这样的：
1.计算输入二维数组的长度
2.按照数组的长度依次选取坐标并与数组中的坐标依次比较，若坐标相同输出一个为零的加速度（质点对自身没有作用力），若坐标不同则给出该坐标所代表的质点所施加的加速度
3.对所有的质点都求过加速度后将加速度求和得到该质点在这一时刻总加速度，（这里体现了mma的一个优点，在mma里，一维数组被视为向量，因此可以直接对分量求和得到总量，相当于Python自带numpy）
4.循环这一过程，得到所有质点的总加速度
5.输出以加速度为元素的数组
这样，我们就把坐标转换成了加速度
1.计算输入二维数组的长度
2.按照数组的长度依次选取坐标并与数组中的坐标依次比较，若坐标相同输出一个为零的加速度（质点对自身没有作用力），若坐标不同则给出该坐标所代表的质点所施加的加速度
3.对所有的质点都求过加速度后将加速度求和得到该质点在这一时刻总加速度，（这里体现了mma的一个优点，在mma里，一维数组被视为向量，因此可以直接对分量求和得到总量，相当于Python自带numpy）
4.循环这一过程，得到所有质点的总加速度
5.输出以加速度为元素的数组
这样，我们就把坐标转换成了加速度
在构造f函数时有三点需要注意，第一是要规定输入的形式，也就是数组的深度和类型，我们的输入是二维实数数组，所以应该规定为{x，_Real,2},其次是要注意统一格式，一开始我规定i==j时If函数返回0，但是0无法与向量求和，因此返回值应该是一维数组{0,0,0}，最后是加速度方向，因为高中时期的固有印象，我没有加负号，结果引力变成了斥力，一画轨发现道体系直接解体
有了加速度，接下来就是解微分方程组了，多体运动的方程组是典型的显示二阶方程组，但是欧拉法只能解一阶的微分方程组，所以我们需要把二阶方程转化为一阶的
这样，我们就得到了：
v（n+1）=v（n）+h*a（n）
x（n+1）=x（n）+h*v（n）
但是a和v在时间间隔h内并非恒定不变的，对于平滑的过程，我们可以给出一个修正：del=（48+20*h+6*h^2+h3）*h/48,所以我们最终得到的结果是：
v（n+1）=v（n）+del*a（n）
x（n+1）=x（n）+del*v（n）
我们在每个时间间隔的起点输入v和x的二维数组，在末尾得到了新的v和x的二维数组，不断地调用输出作为输入，我们就实现了体系各个质点的速度和坐标的等时间间隔离散演化
这样，我们就得到了：
v（n+1）=v（n）+h*a（n）
x（n+1）=x（n）+h*v（n）
但是a和v在时间间隔h内并非恒定不变的，对于平滑的过程，我们可以给出一个修正：del=（48+20*h+6*h^2+h3）*h/48,所以我们最终得到的结果是：
v（n+1）=v（n）+del*a（n）
x（n+1）=x（n）+del*v（n）
我们在每个时间间隔的起点输入v和x的二维数组，在末尾得到了新的v和x的二维数组，不断地调用输出作为输入，我们就实现了体系各个质点的速度和坐标的等时间间隔离散演化
程序在每一步的末尾通过AppendTo函数函数将速度和坐标添加到v和x变量的尾部，每一步的计算结果就以三维数组第一层元素的形式保留了下来。图中共计算了5000步（耗时0.5625秒），v、x的长度就是5001。
计算结束后我们得到的是时间上离散的质点组的运动状态，如果要画天体的运动轨迹，我们需要单个质点在时间上离散的运动状态，也就是说我们需要把质点组分开，在mma里面恰巧有一个函数符合要求，它就是Transpose（矩阵转置），通过它我们可以互换矩阵的行列元素，而x是一个三维数组，我们可以将x视为一个元素为向量的矩阵，这样对x进行转置就相当于把空间上的分组改成了时间上的分组，质点的运动状态也就被分离了出来
计算结束后我们得到的是时间上离散的质点组的运动状态，如果要画天体的运动轨迹，我们需要单个质点在时间上离散的运动状态，也就是说我们需要把质点组分开，在mma里面恰巧有一个函数符合要求，它就是Transpose（矩阵转置），通过它我们可以互换矩阵的行列元素，而x是一个三维数组，我们可以将x视为一个元素为向量的矩阵，这样对x进行转置就相当于把空间上的分组改成了时间上的分组，质点的运动状态也就被分离了出来
有了数据，现在就可以画图了，这也是mma的一大优势，mma拥有众多的数据可视化函数，通过灵活的组合可以画出各种漂亮的图表和动画（画风非常舒服）。这里，我通过Line函数连接单个质点离散的坐标点将其变成轨迹，然后用Graphics3D在三维空间中将轨迹展示出来
我们做的是模拟，既然是模拟，肯定是无法与真实情况完全符合的，必然有误差存在，对于我们目前的程序，这个误差有多大呢？
因为我们模拟的是三体系统，而求解三体系统的解析解是十分困难的，很难直观地判断模拟结果是否精确。所以，为了直接展示程序的精度，我进行了一次两体系统的模拟。我们知道两体系统是稳定且有明确的解析解的，在经典力学的框架下，一个稳定的两体系统的轨道必然是固定的两个椭圆轨道。那么，模拟得到的轨道是怎样的呢
PS：在这次的模拟中，时间步长取0.01，共模拟了1000步
因为我们模拟的是三体系统，而求解三体系统的解析解是十分困难的，很难直观地判断模拟结果是否精确。所以，为了直接展示程序的精度，我进行了一次两体系统的模拟。我们知道两体系统是稳定且有明确的解析解的，在经典力学的框架下，一个稳定的两体系统的轨道必然是固定的两个椭圆轨道。那么，模拟得到的轨道是怎样的呢
PS：在这次的模拟中，时间步长取0.01，共模拟了1000步