The N-Body Problem
可以看到，轨道基本上是稳定的，但是似乎无法完全闭合，轨道的终点稍稍偏离了之前的路径，如果我们延长模拟时间，比如5000步，那么问题将更加明显
为了更清晰地体现这种误差，我将这个两体系统的其中一个星球的质量增加到了原来的十倍，相对速度增加到了原来的2.75倍，同样是5000步
这是50000步的轨道
很明显，星球在不断地沿顺时针“进动”（和相对论没有任何关系）
这是50000步的轨道
很明显，星球在不断地沿顺时针“进动”（和相对论没有任何关系）
误差从何而来呢？
数据是准确的，mma的计算精度是机器精度，f函数严格服从万有引力定律，那么剩下的就只有欧拉法了。作为一种数值求解方法，欧拉法带有二阶的截断误差，它的误差与求解时所采用的步长有关，为O(h^2)，也就是说，每一步的误差为一个很小的数乘上步长的平方。对于0.01的步长来说，O(h^2)是很小的，但是多体系统是一种十分敏感的系统，随着每一步的误差累计，最终的误差就十分明显了。
数据是准确的，mma的计算精度是机器精度，f函数严格服从万有引力定律，那么剩下的就只有欧拉法了。作为一种数值求解方法，欧拉法带有二阶的截断误差，它的误差与求解时所采用的步长有关，为O(h^2)，也就是说，每一步的误差为一个很小的数乘上步长的平方。对于0.01的步长来说，O(h^2)是很小的，但是多体系统是一种十分敏感的系统，随着每一步的误差累计，最终的误差就十分明显了。
看样子，只要减小时间步长就能极大地消除误差了，如果我们把步长设成0.001，那么每一步的误差就会减小一百倍，虽然所需的步数增加了，但最终还是能获得十倍的全局精度。
这是步长为0.001，一万步的轨道
相比于0.01步长一千步的轨道，精度有了极大的提升。但是这种方法是有局限的，虽然我们可以不断地减小时间步长以取得更高的精度，但这是以牺牲计算时间来换取的，在数值模拟中，每一步所需的时间基本上是相同的，更小的时间步长意味着在模拟时长相同的情况下所用的步数更多，所需的计算资源也就增加了。
随着如今电费和CPU越来越贵，这种以时间等比例换取精度的做法无疑会被你的老板和同学打死
这是步长为0.001，一万步的轨道
相比于0.01步长一千步的轨道，精度有了极大的提升。但是这种方法是有局限的，虽然我们可以不断地减小时间步长以取得更高的精度，但这是以牺牲计算时间来换取的，在数值模拟中，每一步所需的时间基本上是相同的，更小的时间步长意味着在模拟时长相同的情况下所用的步数更多，所需的计算资源也就增加了。
随着如今电费和CPU越来越贵，这种以时间等比例换取精度的做法无疑会被你的老板和同学打死
同时，在模拟过程中我还发现了另一个问题，计算时间随计算步数的指数式增加
理论上，欧拉法每一步所需的计算量是相同的，程序的时间复杂度为O(t)，消耗的时间也应该是均匀的，但是很明显，程序所消耗的时间随着步数的增加在快速上升
理论上，欧拉法每一步所需的计算量是相同的，程序的时间复杂度为O(t)，消耗的时间也应该是均匀的，但是很明显，程序所消耗的时间随着步数的增加在快速上升
要获得长时间的准确的轨道数据，这两个问题都是必须要解决的。
首先是误差问题，我们已经知道误差的根本来源是欧拉法，并且单纯的减小时间步长是行不通的，所以只能从算法本身着手。
欧拉法属于龙格库塔法(Runge-Kutta methods)，是一种一阶的也是唯一自洽的一阶显示龙格库塔法。虽然欧拉法运用简单且易于理解，但在实际运用中误差往往难以控制，目前运用最广泛的是经典四阶龙格库塔法(RK4)
其中：
k1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；
k3也是中点的斜率，但是这次采用斜率k2决定y值；
k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值。
龙格库塔法的系数由布切尔表(Butcher tableau)决定，布切尔表是基于泰勒展开推导得到的，通过解一组线性方程使方程的系数与泰勒展式的系数相等，就可以得到一组特解，这组特解对应的就是布切尔表的系数。
对于一阶问题，特解只有一个，就是欧拉法，而对于高阶问题，特解不止一个，因此对应的衍生方法也就多种多样，对于经典RK4，它所对应的特解是最简单的，因此运用也最广泛。
RK4法是四阶方法，也就是说每步的误差是5阶，而总积累误差为4阶。
首先是误差问题，我们已经知道误差的根本来源是欧拉法，并且单纯的减小时间步长是行不通的，所以只能从算法本身着手。
欧拉法属于龙格库塔法(Runge-Kutta methods)，是一种一阶的也是唯一自洽的一阶显示龙格库塔法。虽然欧拉法运用简单且易于理解，但在实际运用中误差往往难以控制，目前运用最广泛的是经典四阶龙格库塔法(RK4)
其中：
k1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；
k3也是中点的斜率，但是这次采用斜率k2决定y值；
k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值。
龙格库塔法的系数由布切尔表(Butcher tableau)决定，布切尔表是基于泰勒展开推导得到的，通过解一组线性方程使方程的系数与泰勒展式的系数相等，就可以得到一组特解，这组特解对应的就是布切尔表的系数。
对于一阶问题，特解只有一个，就是欧拉法，而对于高阶问题，特解不止一个，因此对应的衍生方法也就多种多样，对于经典RK4，它所对应的特解是最简单的，因此运用也最广泛。
RK4法是四阶方法，也就是说每步的误差是5阶，而总积累误差为4阶。
和欧拉法一样，RK4也只能解一阶微分方程，所以我们同样需要降阶，把原方程分解为两个一阶方程并交叉求解，具体步骤如图：
采用RK4后，取0.01的时间步长，5000步，轨迹如下：
可以看到，非常稳，比0.001步长的欧拉法还稳，简直是稳如老狗，那么继续增加步数会如何呢
采用RK4后，取0.01的时间步长，5000步，轨迹如下：
可以看到，非常稳，比0.001步长的欧拉法还稳，简直是稳如老狗
，那么继续增加步数会如何呢
想要继续增加步数，就要解决随步数增加计算时间指数式上升的问题。
这个问题的根源出在变量的大小上，为了保存计算数据，我使用了AppendTo函数，这个函数的作用是将一个元素放在某一数组的尾部，所以每一次保存数据时都需要调用当前步的数据xpres以及之前的所有数据x，对于两体体系，xpres只有48B，但x却会随着计算而累积，对于5000步的计算，x的大小会超过200kB，x越大，调用x所需的时间也就越多，当x继续增加时，CPU的缓存最终将无法完全存储x，这时x将被存储在内存上，调用所需的时间也就进一步延长了
这个问题的根源出在变量的大小上，为了保存计算数据，我使用了AppendTo函数，这个函数的作用是将一个元素放在某一数组的尾部，所以每一次保存数据时都需要调用当前步的数据xpres以及之前的所有数据x，对于两体体系，xpres只有48B，但x却会随着计算而累积，对于5000步的计算，x的大小会超过200kB，x越大，调用x所需的时间也就越多，当x继续增加时，CPU的缓存最终将无法完全存储x，这时x将被存储在内存上，调用所需的时间也就进一步延长了
为了解决这个问题，我采用了截断变量的方法，也就是在变量增长到足够大之前将变量输出为文件然后重置变量。
这种方法的优点是能够保存数据，特别是对于长时间的计算，即使计算被中断，之前已输出的数据是不会丢失的，而且可以通过生成的文件大致估计算了多少。当然缺点也是显而易见的，每一次输出文件都需要时间，这会降低程序的效率，但是相对于之前的情况，这显然利大于弊。
PS1：后来我发现我简直**，用Join函数就能够避开对x函数的调用，不过变量过大时仍然会有降速现象，所以这种分步输出的方式我还是保留了下来。
PS2:在使用Join函数后，RK4（计算量为欧拉法的4倍），50000步所需的时间减少为28s，是使用AppendTo函数时的十分之一
这种方法的优点是能够保存数据，特别是对于长时间的计算，即使计算被中断，之前已输出的数据是不会丢失的，而且可以通过生成的文件大致估计算了多少。当然缺点也是显而易见的，每一次输出文件都需要时间，这会降低程序的效率，但是相对于之前的情况，这显然利大于弊。
PS1：后来我发现我简直**，用Join函数就能够避开对x函数的调用，不过变量过大时仍然会有降速现象，所以这种分步输出的方式我还是保留了下来。
PS2:在使用Join函数后，RK4（计算量为欧拉法的4倍），50000步所需的时间减少为28s，是使用AppendTo函数时的十分之一
虽然轨道看似已经稳定了，但是如果仔细观察，你会发现50000步的轨道依然有微小的偏移，这些偏移自然是来源于计算的误差。为了准确地体现误差的大小，我选取了体系总能量作为参考量。
为了计算体系总能量，我们首先要计算体系的势能（由坐标与质量决定）和动能（由速度与质量决定），分别由fEp与fEk函数求出：
分别把x和v映射到这两个函数中，得到体系每一步的势能和动能的数组，再对应相加就得到了体系总能量的离散数据
用ListPlot函数绘制图像，这是步长0.01，5800步（10周期）的体系能量变化趋势，体系总能量呈阶梯式下降
放大第一周期区间，可以看出能量主要在第290步附近剧烈变化，这段对应的时间正好是两颗星球距离最近的时刻，或者说是星球的速度最快，加速度最大的时刻。
继续放大变化区间到290±20的范围，能量先升后降的趋势更加明显，而且下降区间比上升区间更长更陡，这就导致了每一次经过轨道近点，体系总能量就有一1.2e-7左右的损失。
为了计算体系总能量，我们首先要计算体系的势能（由坐标与质量决定）和动能（由速度与质量决定），分别由fEp与fEk函数求出：
分别把x和v映射到这两个函数中，得到体系每一步的势能和动能的数组，再对应相加就得到了体系总能量的离散数据
用ListPlot函数绘制图像，这是步长0.01，5800步（10周期）的体系能量变化趋势，体系总能量呈阶梯式下降
放大第一周期区间，可以看出能量主要在第290步附近剧烈变化，这段对应的时间正好是两颗星球距离最近的时刻，或者说是星球的速度最快，加速度最大的时刻。
继续放大变化区间到290±20的范围，能量先升后降的趋势更加明显，而且下降区间比上升区间更长更陡，这就导致了每一次经过轨道近点，体系总能量就有一1.2e-7左右的损失。
这是步长0.02，50000步的轨道，轨道远点下降明显
这是因为能量损失主要发生在近点，因为近点的坐标相对固定，所以能量损失主要体现在动能上，导致星球在近点的速度减小，轨道远点随之下降。
这样的体系显然是不满足能量守恒关系的，为了保持体系总能量的稳定，我加入了能量守恒约束函数与动量守恒约束函数：
在每一步运算完成后，先计算体系总能量，与初始能量对比，然后通过调整动能使体系总能量守恒，体系能量守恒之后计算体系动量，与初始动量比较得到一个差值，将这个差值作用到每一个星球的速度上（相当于给参考系一个很小的速度增量），使体系的动量与初值一致
步长0.02,50000步，开启约束的轨迹
这是因为能量损失主要发生在近点，因为近点的坐标相对固定，所以能量损失主要体现在动能上，导致星球在近点的速度减小，轨道远点随之下降。
这样的体系显然是不满足能量守恒关系的，为了保持体系总能量的稳定，我加入了能量守恒约束函数与动量守恒约束函数：
在每一步运算完成后，先计算体系总能量，与初始能量对比，然后通过调整动能使体系总能量守恒，体系能量守恒之后计算体系动量，与初始动量比较得到一个差值，将这个差值作用到每一个星球的速度上（相当于给参考系一个很小的速度增量），使体系的动量与初值一致
步长0.02,50000步，开启约束的轨迹
对于两体系统，要解决能量不守恒问题非常简单，只需要减小步长即可，因为RK4的误差是四阶的，所以步长每减小十倍，单步误差就会减小十万倍，总体误差则减小一万倍。对于目前的例子，步长为0.001时体系的能量波动就会降低到机器精度之下，也就是说已经达到了数值方法的理论极限，无法继续优化了，在图像上体现为能量不变，轨道永久稳定。
但是在多体系统中，天体的轨迹是混沌的，一旦出现近距离的引力弹弓现象，参与该过程的天体的能量波动就会急剧增加，迅速打破体系的能量守恒，严重时甚至会把体系总能量强行抬过零点（多体系统能够稳定存在的必要条件是体系总能量小于零），导致体系直接解体。这就是为什么很多人想要利用一些天文方面的软件和游戏（比如宇宙沙盘）模拟三体系统却往往失败的原因。
多体引力系统的这种能量敏感性来源于引力的物理形式，对于传统的数值模拟方法，比如经典分子动力学模拟，原子的势函数形式使得原子在受力较大时速度较小，速度较大时受力较小，这就保证了运动方程的平滑性（我们知道龙格库塔法基于泰勒展开，而泰勒展开对函数的连续性有要求，因此龙格库塔法在处理平滑问题时非常高效）。而引力则相反，在引力系统中，天体受力大时速度也大，因此在轨道近点运动方程不平滑，这就是误差的最主要来源。好在RK4能够通过减小步长快速提高平滑性，因此我们可以用较小的计算代价来稳定轨道近点处的能量。
要做到这一点，我们需要可变步长法。
多体引力系统的这种能量敏感性来源于引力的物理形式，对于传统的数值模拟方法，比如经典分子动力学模拟，原子的势函数形式使得原子在受力较大时速度较小，速度较大时受力较小，这就保证了运动方程的平滑性（我们知道龙格库塔法基于泰勒展开，而泰勒展开对函数的连续性有要求，因此龙格库塔法在处理平滑问题时非常高效）。而引力则相反，在引力系统中，天体受力大时速度也大，因此在轨道近点运动方程不平滑，这就是误差的最主要来源。好在RK4能够通过减小步长快速提高平滑性，因此我们可以用较小的计算代价来稳定轨道近点处的能量。
要做到这一点，我们需要可变步长法。
更新内容有以下几项
1.加入了脚本输入功能，现在程序可以从指定的文件夹读取质量、尺寸、初始坐标、初始速度数据，避免了更改模拟对象时的麻烦
2.将模拟对象分成了引力体和非引力体两组，非引力体不对其他天体施加引力，对于存在n个引力体，m个非引力体的系统，之前的代码计算量是（n+m）*（n+m-1），现在降低为n*（n+m-1），减小了计算量，尤其是在模拟小行星群时计算速度会有显著的提升
3.引入真实单位，通过设置空间缩放系数与时间缩放系数，将原本非常大的天文单位降低到正常范围（大约在0.1~1000），避免了运动方程的刚性带来的误差
4.引入可变步长法，提升了总体的速度与精度
1.加入了脚本输入功能，现在程序可以从指定的文件夹读取质量、尺寸、初始坐标、初始速度数据，避免了更改模拟对象时的麻烦
2.将模拟对象分成了引力体和非引力体两组，非引力体不对其他天体施加引力，对于存在n个引力体，m个非引力体的系统，之前的代码计算量是（n+m）*（n+m-1），现在降低为n*（n+m-1），减小了计算量，尤其是在模拟小行星群时计算速度会有显著的提升
3.引入真实单位，通过设置空间缩放系数与时间缩放系数，将原本非常大的天文单位降低到正常范围（大约在0.1~1000），避免了运动方程的刚性带来的误差
4.引入可变步长法，提升了总体的速度与精度
可变步长法的关键在于步长的预判，这里我以能量为标准来确定，因此在计算参数里添加了能量收敛判据（Energy Convergence Criterion）
这是可变步长部分的代码
因为可变步长需要多次使用RK4法，因此我把RK4过程也预编译了一遍
RK4函数的输入为{坐标，速度，步长}，输出也为{坐标，速度，步长}，因此可以利用函数嵌套的方式代替循环，
这里我使用Nest函数进行函数迭代，相当于RK4（RK4（RK4(……RK4(x,v,h)））
步长的预判步骤如下：
1.由上一步的分割数cut确定本次的步长h/cut，计算一组步长为h/cut，共cut步的运动过程（该过程的总时长自然为h）
2.得到输出后计算体系中每一个天体单位质量的能量
3.然后以相同的输入数据计算一组步长为h/2cut，共2cut步的运动过程
4.计算单位能量并与之前的总能量做差，得到一个能量差列表
5.选取列表中绝对值最大的数作为本次计算的能量波动demax
6.因为RK的精确度为四阶，因此以demax/ECC四次方根作为调整参数adjust
7.向上取整adjust*cut得到下一步的cut值
这是可变步长部分的代码
因为可变步长需要多次使用RK4法，因此我把RK4过程也预编译了一遍
RK4函数的输入为{坐标，速度，步长}，输出也为{坐标，速度，步长}，因此可以利用函数嵌套的方式代替循环，
这里我使用Nest函数进行函数迭代，相当于RK4（RK4（RK4(……RK4(x,v,h)））
步长的预判步骤如下：
1.由上一步的分割数cut确定本次的步长h/cut，计算一组步长为h/cut，共cut步的运动过程（该过程的总时长自然为h）
2.得到输出后计算体系中每一个天体单位质量的能量
3.然后以相同的输入数据计算一组步长为h/2cut，共2cut步的运动过程
4.计算单位能量并与之前的总能量做差，得到一个能量差列表
5.选取列表中绝对值最大的数作为本次计算的能量波动demax
6.因为RK的精确度为四阶，因此以demax/ECC四次方根作为调整参数adjust
7.向上取整adjust*cut得到下一步的cut值
这种方法的优点是收敛速度快且效率较高，通常在两步以内cut值就能够收敛到稳定区间，对于cut步的运算总共需要计算3cut次，有效计算为2cut次。
当然，问题也是有的，因为每一步的cut值来源于上一步的计算，因此对轨道的平滑性仍然有要求，极近距离的引力弹弓效应还是非常棘手，为了解决这一问题，我引入了两个判定
当cut大于1000时，将cut的线性增长改为指数型，避免cut过大
当adjust大于2时，用下一步的cut重新计算当前步，减少当前步的误差
当然，问题也是有的，因为每一步的cut值来源于上一步的计算，因此对轨道的平滑性仍然有要求，极近距离的引力弹弓效应还是非常棘手，为了解决这一问题，我引入了两个判定
当cut大于1000时，将cut的线性增长改为指数型，避免cut过大
当adjust大于2时，用下一步的cut重新计算当前步，减少当前步的误差