hoshimi、Runge-Kutta法・解ける。

おかげさまで解けました。ありがとうございました。
惑星の運動方程式は以下の通り。ただし定数は省略してあります。
\vec{v} = \frac{d\vec{p}}{dt}
\frac{d^2\vec{p}}{dt^2} = -\frac{\vec{p}}{|\vec{p}|^3}
で、Euler法だと
\Delta \vec{p} = \vec{v}\Delta t
\Delta \vec{v} = -\frac{\vec{p}}{|\vec{p}|^3}\Delta t
としますが、Runge-Kutta法では
f(\vec{v}) = \vec{v}\Delta t
g(\vec{p}) = -\frac{\vec{p}}{|\vec{p}|^3}\Delta t
とおいて
\vec{p_1} = f(\vec{v})
\vec{v_1} = g(\vec{p})
\Delta \vec{p}= f(\vec{v}+\frac{1}{2}\vec{v_1})\Delta t
\Delta \vec{v}= g(\vec{p}+\frac{1}{2}\vec{p_1})\Delta t
と近似するのが2次のRunge-Kutta法になるようです。これでは全然説明になっていないので詳しくはRunge-Kutta法を解説しているページをご覧ください。(式だけ見せたかったので。)
4次も似たようにpの計算をvを使って、vの計算をpを使って解くことになります。

しかしPHPはメンバ変数やメンバ関数を呼び出すのにアロー演算子を使わなければいけなくて、読みにくいです。関数のオーバーライドや演算子の定義ができるC++はやはりすばらしいと思います。次の課題からC++で書こうかなぁ。(これまでは逆らって(?)PHPで書いてきた。)
cf)mimeTexってなかなか使えますね。今後使っていこうかな。