N-Body问题的数值模拟
概述
N-Body 问题 是一个非常著名的物理问题。在天文学上，它具体指对于多个可以通过引力（基于经典牛顿力学）相互作用的天体，预测其中任意一个天体的运动（即在任意时刻的位置和速度）。在化学中，则可以是预测多个相互作用的分子或原子的运动。
当 N=2 时的二体问题已经完全解决；但对于 N=3 时的三体问题，除了在一些特殊条件下可以得到解析解外，一般的三体问题依然没有解析解。
N-Body 问题的模拟(Solver) 是指对于一组输入初始状态（各粒子的质量、位置和速度），输出在用户指定时刻的状态。我们先关注牛顿力学下二维空间的 N-Body 问题，然后再扩展到三维空间。
从历史上看，第一个特殊研究的三体问题是太阳、地球和月亮间的运动。从现代意义上说，它可以指经典力学或量子力学中三个粒子运动的任何问题。
参考 N-Body 问题 - wikipedia 。
参考 《Introduction to ParallelProgramming》 by Peter Pacheco.
数学建模
通用公式
假设在时刻 t ，行星 q 的位置为 $ s_q(t) $，行星 k 的位置为 $ s_k(t) $，则根据牛顿定律，行星 k 作用于行星 q 上的引力为：
$$ f_{qk}(t) = -\frac {Gm_qm_k} { { ( s_q(t) - s_k(t) ) } ^ 3 } ( s_q(t) - s_k(t) ) $$
其中， G 为 万有引力常量（$ 6.67310^{-11} m^3/(kgs^2) $），m 为 行星质量，s 和 f 均为带方向的向量。
进一步，我们可以计算所有其他行星作用于行星 q 上的引力和（向量和）：
$$ F_q(t) = \sum_{k=0,k\neq q}^{n-1}f_{qk} = -Gm_q \sum_{k=0,k \neq q} ^ {n-1} \frac {m_k} { { ( s_q(t) - s_k(t) ) } ^ 3 } [s_q(t) - s_k(t)] $$
最后，我们可以计算行星 q 在时刻 t 时的加速度：
$$ a_q(t) = -G \sum_{j=0,j \neq q}^{n-1} \frac {m_j} { { ( s_q(t) - s_j(t) ) } ^ 3 } [s_q(t) - s_j(t)] $$
实际上，也可以用下面的等式表示：
$$ a_q(t) = {s_q}’’(t) = \frac {F_q(t)} {m_q} $$
数值模拟
如果要求 $ s_q(t) $，直接求解上节最后的微分公式是困难的，但是可以使用数值模拟方法来求得其近似值。我们使用最简单的欧拉方法，作为数值模拟方法。对于在定义域内连续的函数 $ g(t) $，在其曲线上有点$ (t_0, g(t_0)) $，那么当 $ t = t_0 + \Delta{t} $时，
$$ y = g(t_0) + g’(t_0)(\Delta{t}) $$
且当 $ \Delta{t} $足够小时：
$$ g(t_0 + \Delta(t)) \approx g(t_0) + g’(t_0)(\Delta{t}) $$
以此，我们可以模拟 $ s_q(t) $和$ v_q(t) $的值。假设，在时刻 0，行星 q 的位置和速度为 $ s_q(0) $ 和 $ v_q(0) $，那么，在时刻 $ \Delta{t} $，其位置和速度为：
$$ s_q(\Delta{t}) \approx s_q(0) + \Delta{t}s_q’(0) = s_q(0) + \Delta{t}v_q(0) $$
$$ v_q(\Delta{t}) \approx v_q(0) + \Delta{t}v_q’(0) = v_q(0) + \Delta{t}a_q(0) = v_q(0) + \Delta{t}\frac{1}{m_q}F_q(0) $$
当我们求得$s_q(\Delta{t})$ 和 $v_q(\Delta{t})$后，我们可以使用相同的公式计算 $s_q(2\Delta{t})$ 和 $v_q(2\Delta{t})$，以此类推。