The N-Body Problem
在目前的基础上，我进行了一次102体的模拟，引力体为地月系，非引力体是以地月拉格朗日点L4为中心，标准差为一万公里，呈正态分布的一百颗小行星，速度同样为正态分布，标准差为0.03km/s（月球公转速度1.023km/s）
本次模拟耗时9816秒，时间步长为1200秒，一共模拟了30个循环，每个循环600步，共18000步，总时长2160万秒（约250天），采样循环为5，采集了3600帧，
ECC为10^-13，小于总能量的百万亿分之一，可以认为模拟是精确的
我记录了每一步的cut值，这是绝大部分cut值的分布
这是所有cut值的分布，其中有部分cut值超出了1000，这意味着这些步的能量误差超过了10^-13
为了方便观察把点连了起来，可以看到最大的步数在12500左右（实际应取的步数约为156000），这些峰使得计算过程中实际的总步数达到了346942，远远超过了基础步数18000
造成这种极大cut值的根本原因在于基本步长的分割方法，我采用的方法是均匀分割法，这就意味着一旦一个基本步中出现了天体间的近距离相互作用，它就需要被分割成能够满足近距离相互作用的步长，而实际上仅有一小部分时间需要如此短的步长，其余时间的大部分步数就被浪费掉了，我们陷入了与固定步长时相似的困境
ECC为10^-13，小于总能量的百万亿分之一，可以认为模拟是精确的
我记录了每一步的cut值，这是绝大部分cut值的分布
这是所有cut值的分布，其中有部分cut值超出了1000，这意味着这些步的能量误差超过了10^-13
为了方便观察把点连了起来，可以看到最大的步数在12500左右（实际应取的步数约为156000），这些峰使得计算过程中实际的总步数达到了346942，远远超过了基础步数18000
造成这种极大cut值的根本原因在于基本步长的分割方法，我采用的方法是均匀分割法，这就意味着一旦一个基本步中出现了天体间的近距离相互作用，它就需要被分割成能够满足近距离相互作用的步长，而实际上仅有一小部分时间需要如此短的步长，其余时间的大部分步数就被浪费掉了，我们陷入了与固定步长时相似的困境
这次更新主要改进了可变步长法，用一个自定义函数split代替了原来通过cut值来分割步长的方式，代码立刻简洁了很多
split函数（上一楼的split有点问题，每步的ECC没变，会导致累计误差变大）
这个函数的基本思想是自适应递归算法，在一个基本步中，函数会将输入的时间步长按照与之前类似的方法处理一遍（先计算一个步长为h的步，然后从相同的起点计算两个步长为h/2的步，分别计算体系总能得到误差），如果能量误差在收敛范围之内，函数会直接返回两步计算的结果（两步计算的误差只有单步的1/16），如果不在收敛范围内，函数会以原步长的二分之一为参数调用自身，得到结果后再以相同的步长调用这个结果，相当于将原先的一步分为了两步，每步所要满足的能量误差也会减半，使每一个基本步的总体能量误差保持在ECC以内。这种“下探”的方法会消耗更多的试探步数（原来的算法有1/3的步数是用于试探的，现在则变为1/3~1/2），但不会有被浪费的步数，在轨道不稳定时更加高效。但实际计算的结果表明两者的运算速度差别不大，这可能是因为这种递归函数需要多次调用自身，拖慢了总体的速度。
PS：之所以使用动态函数而不预编译split是因为mma无法编译带@@（Apply）的函数，深层次的原因我暂时也不清楚
split函数（上一楼的split有点问题，每步的ECC没变，会导致累计误差变大）
这个函数的基本思想是自适应递归算法，在一个基本步中，函数会将输入的时间步长按照与之前类似的方法处理一遍（先计算一个步长为h的步，然后从相同的起点计算两个步长为h/2的步，分别计算体系总能得到误差），如果能量误差在收敛范围之内，函数会直接返回两步计算的结果（两步计算的误差只有单步的1/16），如果不在收敛范围内，函数会以原步长的二分之一为参数调用自身，得到结果后再以相同的步长调用这个结果，相当于将原先的一步分为了两步，每步所要满足的能量误差也会减半，使每一个基本步的总体能量误差保持在ECC以内。这种“下探”的方法会消耗更多的试探步数（原来的算法有1/3的步数是用于试探的，现在则变为1/3~1/2），但不会有被浪费的步数，在轨道不稳定时更加高效。但实际计算的结果表明两者的运算速度差别不大，这可能是因为这种递归函数需要多次调用自身，拖慢了总体的速度。
PS：之所以使用动态函数而不预编译split是因为mma无法编译带@@（Apply）的函数，深层次的原因我暂时也不清楚
本来我这次还有许多改进的思路，比如去掉f函数中的判断、用Sum函数代替Total+Table、将函数变为并行模式等等，但最后发现大部分改进还不如原来的快，而且存在很大的问题，因此只保留了可变步长法的修改。
PS：mma的并行运算是真的辣鸡（也可能是我太菜了）
PS：mma的并行运算是真的辣鸡
（也可能是我太菜了）
之前的split函数有BUG，每一次下探返回的ECC没有乘二，导致ECC只减不增，下探的次数一多就收敛不了，同时为了防止ECC被切分到舍入误差以下我给它设了一个1e-15的基础值
在修改过的split函数基础上我进行了一次长时间的模拟，耗时30996.7秒
仍然是之前的地月系统+L4附近100小行星的体系，共六万步（之前是18000步），从时间上看两者的差别不大。
前3600步计算所需的实际步数：
前18000步：
两者差别不大，相对而言split函数长时间模拟占优势，（实际上均匀分割法应消耗更多的步数，因为当cut超过1000时程序对cut值进行了限制）
从步数分布上来看（蓝色均匀分割，黄色split函数），split函数的波动更小，但在天体的运动更平缓时，split函数并不占优势（需要更多的试探步），且因为需要多次调用动态函数，单步的速度也慢于平均分割法
仍然是之前的地月系统+L4附近100小行星的体系，共六万步（之前是18000步），从时间上看两者的差别不大。
前3600步计算所需的实际步数：
前18000步：
两者差别不大，相对而言split函数长时间模拟占优势，（实际上均匀分割法应消耗更多的步数，因为当cut超过1000时程序对cut值进行了限制）
从步数分布上来看（蓝色均匀分割，黄色split函数），split函数的波动更小，但在天体的运动更平缓时，split函数并不占优势（需要更多的试探步），且因为需要多次调用动态函数，单步的速度也慢于平均分割法
放大前3600步的区间
主要分布范围
可以看到split函数确实更加平稳，均匀分割法会导致cut值上下震荡。
放大1400~1600区间，震荡更加明显
主要分布范围
可以看到split函数确实更加平稳，均匀分割法会导致cut值上下震荡。
放大1400~1600区间，震荡更加明显
总体而言，在相同的ECC下，split函数的精度更高（因为每一步的能量都是严格收敛的），在极端状况下的表现更好，代价是需要更多的试探步数，不适合稳定体系
到目前为止，虽然程序在精度上已经足够让人满意，可以在较短的时间内给出单步能量波动在十亿分之一以内的长时间模拟的结果，但是还存在一个硬伤——没有考虑天体的碰撞。
我从采集到的的12000帧里提取出了所有小行星与地球和月球最近的十个距离（单位为一万公里）
地球：
月球：
月球的半径为1737km，也就是说有1~4颗小行星在模拟过程中与月球的距离曾经小于月球半径，换句话说，与月球发生了撞击。
缺少碰撞判断一是不符合客观现实，其次也会导致计算量增加（天体距离越近相互作用越强越快，误差也就越大）
到目前为止，虽然程序在精度上已经足够让人满意，可以在较短的时间内给出单步能量波动在十亿分之一以内的长时间模拟的结果，但是还存在一个硬伤——没有考虑天体的碰撞。
我从采集到的的12000帧里提取出了所有小行星与地球和月球最近的十个距离（单位为一万公里）
地球：
月球：
月球的半径为1737km，也就是说有1~4颗小行星在模拟过程中与月球的距离曾经小于月球半径，换句话说，与月球发生了撞击。
缺少碰撞判断一是不符合客观现实，其次也会导致计算量增加（天体距离越近相互作用越强越快，误差也就越大）
发现了一个bug，今天上午我试着模拟了一下比邻星（带行星）略过半人马A、B系统的情形，时间步长43200秒，每个循环5000步。刚开始的十几个循环平均15s一个，然后我就把程序挂着去恰饭，结果回来一看还卡在那一步。中断一查，最后一步下探了两百多万次还没结束。
仔细检查了一遍程序，我发现问题出在势能函数上：
求势能需要除以天体之间的距离，距离本身存在截断误差和舍入误差，当距离很小时，势函数具有刚性（一个很小的输入误差会导致很大的输出误差），当天体间的距离足够小且速度较快时，函数的刚性会使能量的总误差随步长减小的速度无法在ECC足够大时追上ECC的减小速度，导致split函数无穷递归，收不回来。
为了解决这个问题，我放弃了将能量作为收敛判据，改为以坐标和速度作为收敛判据（我感觉我之前脑抽了，用能量做收敛判据还需要每次下探求三次总能量，吃力不讨好）。
这是新的split函数：
每一步以坐标和速度误差的模作为判断条件，与XCC和VCC比较。
实际上，以能量为收敛判据也是可以的，只不过需要放弃一部分精度来增加误差的容限，防止出现误差减小追不上收敛标准减小的情况。
仔细检查了一遍程序，我发现问题出在势能函数上：
求势能需要除以天体之间的距离，距离本身存在截断误差和舍入误差，当距离很小时，势函数具有刚性（一个很小的输入误差会导致很大的输出误差），当天体间的距离足够小且速度较快时，函数的刚性会使能量的总误差随步长减小的速度无法在ECC足够大时追上ECC的减小速度，导致split函数无穷递归，收不回来。
为了解决这个问题，我放弃了将能量作为收敛判据，改为以坐标和速度作为收敛判据（我感觉我之前脑抽了，用能量做收敛判据还需要每次下探求三次总能量，吃力不讨好）。
这是新的split函数：
每一步以坐标和速度误差的模作为判断条件，与XCC和VCC比较。
实际上，以能量为收敛判据也是可以的，只不过需要放弃一部分精度来增加误差的容限，防止出现误差减小追不上收敛标准减小的情况。
这次模拟一共有五十万个基础步，实际步数2083972步，每一步的切分次数如下图：
从这张图可以很明显地看出，大部分步数都集中在3,7,15,63等几个特殊值上，这是因为本次设置的四体体系的运动轨迹相对平缓，大部分基本步内的速度变化不大，所以都被等分了，切分次数也就对应于∑2^n，当切分次数超过127时，基本步内的运动不再均匀，分割出来的步长长短不一，因此分割次数大部分都落在了特殊值之间
这是将每一步连接之后的图像，放大之后会发现整个图像由一个个尖锐的峰组成
放大水平区间，每一个峰的基本形状开始显现，从这里可以看出，峰的分布范围很小，大部分情况下基础步是单步的，这也解释了为什么每个基础步的平均切分次数只有4.168
继续放大，可以看到峰的细节，切分次数是非常明显的2指数阶梯式增加
在这一次的模拟中，split函数合理地对步长进行了控制，表现基本符合预期，相较于相同步长的不可变步长情形，它以十二倍的计算量将计算误差降低了大约十六亿倍
从这张图可以很明显地看出，大部分步数都集中在3,7,15,63等几个特殊值上，这是因为本次设置的四体体系的运动轨迹相对平缓，大部分基本步内的速度变化不大，所以都被等分了，切分次数也就对应于∑2^n，当切分次数超过127时，基本步内的运动不再均匀，分割出来的步长长短不一，因此分割次数大部分都落在了特殊值之间
这是将每一步连接之后的图像，放大之后会发现整个图像由一个个尖锐的峰组成
放大水平区间，每一个峰的基本形状开始显现，从这里可以看出，峰的分布范围很小，大部分情况下基础步是单步的，这也解释了为什么每个基础步的平均切分次数只有4.168
继续放大，可以看到峰的细节，切分次数是非常明显的2指数阶梯式增加
在这一次的模拟中，split函数合理地对步长进行了控制，表现基本符合预期，相较于相同步长的不可变步长情形，它以十二倍的计算量将计算误差降低了大约十六亿倍
这是前十五万步的轨道，相当于103年，其中红色的轨迹是半人马A的，蓝色是半人马B，黑色是比邻星，绿色是行星，A、B完全按照真实情形设置（远端距离35.6AU，相对速度4.87km/s，周期79.91年），比邻星现实中距离AB系统约15000AU，为了模拟三体系统我把它移动到了距离AB中心46AU左右的位置，然后在距离比邻星1AU的位置放了一个虚构的行星
这次模拟出现了一个比较有意思的现象，在比邻星第一次掠过半人马A时，行星被半人马A捕获了，然后行星一直围绕半人马A运行，在此后的多次近距离相互作用过程中一直没有被弹出去
这次模拟出现了一个比较有意思的现象，在比邻星第一次掠过半人马A时，行星被半人马A捕获了，然后行星一直围绕半人马A运行，在此后的多次近距离相互作用过程中一直没有被弹出去
花了这么多时间主要是为了保证体系的质量与动量守恒，普通的三体模拟代码（我见过的）会在每一步计算过后执行一遍碰撞判断，一般是扫描一个坐标-坐标的上（下）三角阵，逐次计算天体两两之间的距离，与两个天体的半径之和比较，大于半径之和认为碰撞不发生，小于说明碰撞发生，碰撞发生后生成一个新的天体，执行完所有判断后清除发生过碰撞的天体，或者生成新天体后直接移除参与碰撞的天体（这需要动态更新数组）
通常，这种两做法都是没有问题的，但在极端条件有bug
考虑这种情况：两个天体同时与另一个天体相撞，这时，第一种方法会使处于判断逻辑序列中间的天体被计算两次，使体系凭空增加一个中间天体的质量，第二种方法则有可能使第三个天体避免碰撞，与实际情况不符
通常，这种两做法都是没有问题的，但在极端条件有bug
考虑这种情况：两个天体同时与另一个天体相撞，这时，第一种方法会使处于判断逻辑序列中间的天体被计算两次，使体系凭空增加一个中间天体的质量，第二种方法则有可能使第三个天体避免碰撞，与实际情况不符
为了避免上述问题，我加入了一个分组函数：
这个函数的作用是将所有相连的天体分到同一组，以分组为依据执行碰撞，比如有1,2,3,4,5,6六个天体，其中1,2相连，2,3相连（但1,3不相连），4孤立，5,6相连，那么函数就会将其分为{1,2,3}，{4}，{5,6}三组
这个函数的逻辑过程是这样的：
1.将输入的数列设为原始组
2.将原始组的第一个元素扔到一个新的组里面
3.用原始组的每一个元素与新组的第一个元素匹配并判断，符合要求（发生碰撞）的扔到新组
4.用原始组的每一个元素与新组的下一个元素匹配并判断，符合要求（发生碰撞）的扔到新组
5.重复4，直到新组的最后一个元素
6.重复2，直到原始组没有元素
PS：我不是CS专业的，不清楚这种算法具体叫啥，有懂的朋友麻烦科普一下
这个函数的作用是将所有相连的天体分到同一组，以分组为依据执行碰撞，比如有1,2,3,4,5,6六个天体，其中1,2相连，2,3相连（但1,3不相连），4孤立，5,6相连，那么函数就会将其分为{1,2,3}，{4}，{5,6}三组
这个函数的逻辑过程是这样的：
1.将输入的数列设为原始组
2.将原始组的第一个元素扔到一个新的组里面
3.用原始组的每一个元素与新组的第一个元素匹配并判断，符合要求（发生碰撞）的扔到新组
4.用原始组的每一个元素与新组的下一个元素匹配并判断，符合要求（发生碰撞）的扔到新组
5.重复4，直到新组的最后一个元素
6.重复2，直到原始组没有元素
PS：我不是CS专业的，不清楚这种算法具体叫啥，有懂的朋友麻烦科普一下
为了节约计算资源，我在group函数前面设置了一个判断函数
这个函数耗时约为accelerate函数的70%，它的作用是判断当前帧是否有碰撞发生，如果发生了，将执行group函数并执行碰撞过程
这个函数耗时约为accelerate函数的70%，它的作用是判断当前帧是否有碰撞发生，如果发生了，将执行group函数并执行碰撞过程
这几天debug最主要的对象是split函数：
为了加入碰撞判断，我先是把程序整体的参数传递过程优化了一遍，然后将碰撞判断和碰撞执行耦合进split函数中。
这里最主要的难点是碰撞执行的参数传递过程，因为涉及到递归函数，所以一不小心就会陷入死循环或者参数长度和类型不匹配，再加上mma要debug只能查堆栈（mma确实不适合大规模编程），所以花费了大量时间
为了加入碰撞判断，我先是把程序整体的参数传递过程优化了一遍，然后将碰撞判断和碰撞执行耦合进split函数中。
这里最主要的难点是碰撞执行的参数传递过程，因为涉及到递归函数，所以一不小心就会陷入死循环或者参数长度和类型不匹配，再加上mma要debug只能查堆栈（mma确实不适合大规模编程），所以花费了大量时间
因为数组纬度不是统一的了，所以轨迹处理也更麻烦了，需要通过天体的ID逐个匹配
动画的生成也更加繁琐，尤其是轨迹的绘制需要单独的判断逻辑
动画的生成也更加繁琐，尤其是轨迹的绘制需要单独的判断逻辑
发现了一些比较有趣的现象
我模拟了一个带大卫星的双行星系统
a星1.1地球质量，b星1.4地球质量，a、b距离57306.2km，相对速度4167.42km/s，在孤立情况下恰好能够以稳定的正圆轨道和一天的周期互相绕行，再距离a、b质心38.22万公里处放了一颗1/64.38地球质量，速度1613.7m/s的卫星
这是前五千步a、b的轨道（每步一分钟）：
这是两万五千步的：
我模拟了一个带大卫星的双行星系统
a星1.1地球质量，b星1.4地球质量，a、b距离57306.2km，相对速度4167.42km/s，在孤立情况下恰好能够以稳定的正圆轨道和一天的周期互相绕行，再距离a、b质心38.22万公里处放了一颗1/64.38地球质量，速度1613.7m/s的卫星
这是前五千步a、b的轨道（每步一分钟）：
这是两万五千步的：
仅仅这两张图看不出什么，我又画了a、b距离随时间变化的图像
这是25000步的数据，大约有10km的震荡
这是50000步，可以看到一个清晰地正弦变化趋势
250000步：
2500000步（相当于42.78年），在总体上体现出了一个周期在百年量级的超低频震荡
我本来以为这个是由卫星的潮汐力引起的，结果一画月球的能量变化发现周期对不上：
月球能量变化的周期约为17.36天，与卫星的公转周期吻合，但是a、b距离的二级振荡周期大约只有此周期的一半
这是25000步的数据，大约有10km的震荡
这是50000步，可以看到一个清晰地正弦变化趋势
250000步：
2500000步（相当于42.78年），在总体上体现出了一个周期在百年量级的超低频震荡
我本来以为这个是由卫星的潮汐力引起的，结果一画月球的能量变化发现周期对不上：
月球能量变化的周期约为17.36天，与卫星的公转周期吻合，但是a、b距离的二级振荡周期大约只有此周期的一半
还有就是这颗卫星的能量实际上是不稳定的，虽然能量的均值为-0.00209，在零以下，但是能量的振荡使卫星的总能量有接近一半的时间在零以上