数值问题高精度求解

Ickiverar

数学家们已经从理论上证明，对于特定的初始位置与速度，存在这样一条简单的8字型轨道，三个遵循牛顿引力互相吸引的等质量质点可以在其上进行严格周期的运动。它们的轨迹严格重合，每绕出一个8字即为一个周期，如下图所示：

为了简化问题，假设三个质点的质量都是1，万有引力常量也是1，则引力加速度就是 a = 1/r^2 ，r是两个质点的距离。
考虑这个周期轨道的初值，如下所示：

初值写出来是
r1 = (  1 , 0 ), v1 = (   vx ,   vy );
r2 = (  0 , 0 ), v2 = ( -2vx , -2vy );
r3 = ( -1 , 0 ), v3 = (   vx ,   vy );
问题是，求出对应于图中的严格周期轨道的 vx , vy 值，尽可能的高精度。
我是早就了解了这个问题，从2017年2月2日忽然想要计算这初值的高精度结果，于是先去网上找了vx与vy的近似值
vx~0.347111, vy~0.532728
然后开始算，算了一下午大概得到了40位的精度（不一定对），其中前24位（我认为）可以确定是正确的。
vx~0.3471168881189269382427769203203008662474
vy~0.5327249453880302292620279876919262703549
我现在只有非常弱的算法，每往后算一定位数就要花数倍时间，这是几何式的；
而且虽然用了14阶Runge-Kutta这样高阶的算法，但是系数表却只有网上下载的60位，我不会算它的系数表，这限制了进一步的计算。
一切数学软件都可以用，这是一个实际问题。我用的是Mathematica。

补充内容 (2018-12-15 02:28):
图挂了，8楼有最新的结果和图。

happysxyf

这个只要能解出受力方程的x,y表达式，然后剩下的工作就只是用级数高精度去逼近x,y。我需要你的力学方程，也就是你用的x，y的方程式什么，最后的解的表达式是什么，然后才能分析出可行的高精度算法。
虽然这只是个等质量的三体问题，但是我没研究过，不会列物理方程式。

Ickiverar

@happysxyf
以下的$r$都是矢量。
比如三个球的位置是$r_1,r_2,r_3$，那么：
球2和球3对$r_1$的加速度（对时间的二阶导数）的贡献是线性相加的。
球2对$r_1$的加速度，是沿$r_12$的，大小为$r_12^-2$的矢量，写出来就是${r_12}/{|r_12|}·|r_12|^-2=r_12/{|r_12|^3}$，其中$r_12=r_2-r_1$，球3的贡献类似。
那就有以下关系：
${d^2r_1}/{dt^2}=r_12/{|r_12|^3}+r_13/{|r_13|^3}$
${d^2r_2}/{dt^2}=r_23/{|r_23|^3}+r_21/{|r_21|^3}$
${d^2r_3}/{dt^2}=r_31/{|r_31|^3}+r_32/{|r_32|^3}$
其中$r_{ij}=r_j-r_i$。
初始条件就是
$r_1=(1,0), v_1={dr_1}/{dt}=(v_x,v_y)$
$r_2=(0,0), v_2={dr_2}/{dt}=(-2v_x,-2v_y)$
$r_3=(-1,0), v_3={dr_3}/{dt}=(v_x,v_y)$

happysxyf

Ickiverar 发表于 2017-2-4 13:04
@happysxyf
以下的$r$都是矢量。
比如三个球的位置是$r_1,r_2,r_3$，那么：

这个微分方程你用mathmaticas解的？能不能推导出解的表达式？

Ickiverar

happysxyf 发表于 2017-2-4 13:15
这个微分方程你用mathmaticas解的？能不能推导出解的表达式？

显然不能。
三体问题的解无法写成解析式，甚至无法写成一致收敛级数的和。似乎是庞加莱研究过的命题。

Ickiverar

经过了三天的计算，我可以确信我得到了64位的可信精度：
vx~0.3471168881189269382427769203203008662474073371707833534398872083
vy~0.5327249453880302292620279876919262703549110901882196558216145751
最后一位可能因四舍五入而差1.
用我现在的算法，进一步的计算是不可能的。

happysxyf

Ickiverar 发表于 2017-2-5 09:54
经过了三天的计算，我可以确信我得到了64位的可信精度：
vx~0.34711688811892693824277692032030086624740 ...

好吧，看来这是个棘手的问题，需要有力的数学工具。

Ickiverar

最近想起了这个事情，又捡起来重新算了算。
还是用的12阶RungeKutta，但是我把它的系数表算到了1000位，又搞了许多常数优化方法，比当年的算法和实现加速了大概百分之两三千吧。
另外，新的计算所用的姿态和1楼的设定不一样，1楼的设定是所用初值最少的设定，只用到了两个参数就能搞出来周期轨道，但是轨道算完是斜的......
新的计算考虑了该周期轨道中所有重要的8个参数，包括那个倾斜角 i=0.2455475958...，揭露了轨道的所有对称性（至少有关于x轴和y轴的对称性）。
最终将 8 个重要参数都算到了100位以上，计算用时约16小时：


上图主要给出了两个重要姿态的各球位置和速度。这两个姿态的时间间隔为 T ，而整个系统的周期为 12 T.
从这两个姿态证明轨道的周期性是显然的。
————————————————
另附：在上述计算之前，我也用一楼的姿态设定算了一下：
精度稍高些，有105位（注：这一计算进行的时候，还有一些常数优化没有做，所以耗时60小时。算完才想我为什么不把这个轨道其他参数一块算了呢 orz ）：
vx~0.347116888118926938242776920320300866247407337170783353439887208345636425244992106627630873723412431953247
vy~0.532724945388030229262027987691926270354911090188219655821614575145018704237828103134902339396912268179360
4T~2.10863799430873722586300011132234902109754996996471569739141396536423423362192613713494903683486250883026

wayne

如果我们换种思路： 求这样的轨道曲线方程，使得其上有三点，随着各自运动的时候，始终满足 引力方程。 这个就有点类似于 变分方程了。

Ickiverar

@mathe
3#给出的条件其实就是8#中 t=0 的姿态，只不过旋转了一个角度，使得球（而非轨道长轴）呈水平，这样可以减少一个参数。
从 8# 的两个姿态可以证明周期性成立（但稳定性不能证明。我也不知道如何证明稳定性），故而将初值这样参数化。
————————————————————————
对周期性的证明：

0.设三个球分别为红绿蓝，位置分别为 x1=(x1x,x1y), x2=(x2x,x2y), x3=(x3x,x3y), 速度分别为 v1，v2，v3，均是时间的函数。

1.考虑 8# 中 t=0 的状态，系统重心（三个球位置的平均）为(0,0)，重心速度（三个球速度的平均）为(0,0)；姿态可以由三个参数描述，vx，vy，i，故自由度为3.

2.如 8# 中右图，如果规定 T 为当红球运动到x轴上的时刻（即 x1y=0 ），且要求：
   a. 红球水平速度 v1x=0
   b. 绿蓝两球水平位置相等 x2x=x3x
   c. 绿蓝两球竖直速度相等 v2y=v3y
  那么，根据动量守恒定律，系统的重心和重心速度仍然应是(0,0)，所以一定有:
    x1x=-x2x/2=-x3x/2 （设为A）
    x2y=-x3y  （设为B）
    v1y=-v2y=-v3y   （设为uy）
    v2x=-v3x   （设为ux）
  也就是8#中 t=T 的状态。

注意到 t=0 时有三个参数自由度，而 t=T 时有三个限制，所以恰好可以形成一个适定方程组，解出初值参数 vx，vy，i，使得存在这样一条轨道，从 t=0 时左图的状态演化到 t=T 时右图的状态。
（也可以认为是四个参数自由度（包括T）和四个限制（ x1y=0 ），同样也是适定的）

3. 现在 0~T 的轨道已经证明成立了，剩下的运动可以由对称性得到：
T~2T 的运动，恰好是 0~T 的 时间反演+沿x轴对称+交换蓝绿球
2T~4T 的运动，恰好是 0~2T 的 沿y轴对称+交换三个球
4T~8T 的运动，完全是 0~4T 交换三个球的结果
8T~12T 的运动，同样是 0~4T 交换三个球的结果
t=12T 时，整个系统已经回到 t=0 的状态，故周期性成立。