龙格-库塔（RK4）方法

对于从 $t 到 $t + h 的仿真步骤，我们需要以下输入

• $\vec{x}^t$: 时间 $t 时所有粒子的位置。
• $\vec{v}^t$: 时间 $t 时所有粒子的速度。
• $\vec{a}(\vec{x}^t, \vec{v})$: 一个函数，给定所有位置和速度，计算加速度。

输出，我们将计算 $\vec{x}^{t + h}$ 和 $\vec{v}^{t + h}$ 的近似值。

RK4积分方法的工作原理如下

• 计算 $\vec{a}^t = \vec{a}(\vec{x}^t, \vec{v})$
• 计算 $\vec{x}^{mid_1} = \vec{x}^t + \frac{h}{2}\vec{v}^t$ 和 $\vec{v}^{mid_1} = \vec{v}^t + \frac{h}{2}\vec{a}^t$
• 计算 $\vec{a}^{mid_1} = \vec{a}(\vec{x}^{mid_1}, \vec{v}^{mid_1}))$
• 计算 $\vec{x}^{\text{mid}_2} = \vec{x}^t + \frac{h}{2}\vec{v}^{\text{mid}_1}$ 和 $\vec{v}^{\text{mid}_2} = \vec{v}^t + \frac{h}{2}\vec{a}^{\text{mid}_1}$
• 计算 $\vec{a}^{\text{mid}_2} = \vec{a}((\vec{x}^{\text{mid}_2}, \vec{v}^{\text{mid}_2}))$
• 计算 $\vec{x}^{\text{end}} = \vec{x}^t + h\vec{v}^{\text{mid}_2}$ 和 $\vec{v}^{\text{end}} = \vec{v}^t + h\vec{a}^{\text{mid}_2}$
• 计算 $\vec{a}^{end} = \vec{a}(\vec{x}^{end}, \vec{v}^{end})$
• 近似（离散化误差 $O(h^5)$）
  • $\vec{x}^{t + h} \approx \vec{x}^t + h \frac{\vec{v}^t + 2 \vec{v}^{mid_1} + 2 \vec{v}^{mid_2} + \vec{v}^{end}}{6}$
  • $\vec{v}^{t + h} \approx \vec{v}^t + h \frac{\vec{a}^t + 2 \vec{a}^{mid_1} + 2 \vec{a}^{mid_2} + \vec{a}^{end}}{6}$