導入
ここではケプラー問題
$$
\frac{d^2}{dt^2}{\boldsymbol r} = - \frac{\mu}{r^3} \boldsymbol r
$$
に対して時間変換を施すことで，離心率が髙い場合などでも効率よく数値計算できるようにする方法を紹介していきます．
表記のルールですが，$\boldsymbol r$のような太字はベクトル量，同じ文字が太字ではなくなったときは$ r = | \boldsymbol r |$のようにそのベクトルの大きさを示すものとします．
実時間$t$に対して，仮想時間$s$を
$$
dt = r\, ds
$$
となるように導入します．$s$が等幅で刻まれても，彗星が太陽に近づいたら$t$は細かく刻まれるというイメージです．こうやって$t$についての常微分方程式を$s$についての常微分方程式に移し替えていきます．
時間微分の表記についていくつかの約束事を導入します．
実時間で微分するときは，上にドット（・）をつけます：
$$
\frac{d}{dt}[\, ] = \dot{[\,]}
$$
仮想時間で微分するときは右にプライム（'）をつけます：
$$
\frac{d}{ds}[\, ] = [\,]'
$$
両方の微分の関係は
$$
\frac{d}{ds} = r \frac{d}{dt}
$$
ですから
$$
[\,]' = r \dot{[\,]}
$$
と書くことができます．
$$
t' = r
$$
や
$$
\boldsymbol r' = r\, \dot{\boldsymbol r}
$$
などの書き方ができます．
運動方程式
$\boldsymbol r$の2階微分$\boldsymbol r''$の具体的なかたちを求めていきましょう．仮想時間での2階微分は，
$$
[\,]'' = r \frac{d}{dt} \left( r\, \dot{[\,]} \right) = r^2 \ddot {[\,]} + r \dot r \dot{[\,]}
$$
で計算できます．
$\dot r$の具体的なかたちですが，$\boldsymbol v = \dot {\boldsymbol r}$として速度を用いると
$$
\dot r = \frac{d}{dt} \left( \boldsymbol r \cdot \boldsymbol r \right)^{1/2} = \frac{(\boldsymbol r \cdot \boldsymbol v)}{r}
$$
のようにして計算できます．
これで運動方程式
$$
\ddot {\boldsymbol r} = - \frac{\mu}{r^3} \boldsymbol r
$$
を変換する準備が整いました．変換後の運動方程式は，
$$
\boldsymbol r '' = -\frac{\mu}{r} \boldsymbol r + (\boldsymbol r \cdot \boldsymbol v) \boldsymbol v
$$
になります．これ，実は調和振動子になっています（ただし中心は原点からずれている）．そんなこと言われてもとても信じられませんよね（この式は調和振動子にはとても見えない）．
とりあえず天下り的に離心率ベクトルの定義を貰ってきます．
$$
\boldsymbol e = \left( \frac{v^2}{\mu} - \frac{1}{r} \right) \boldsymbol r - \frac{(\boldsymbol r \cdot \boldsymbol v)}{\mu} \boldsymbol v
$$
ここで証明はしませんが，これが保存量になっています．これを変形すると
$$
(\boldsymbol r \cdot \boldsymbol v) \boldsymbol v = \left( v^2 - \frac{\mu}{r} \right) \boldsymbol r - \mu \boldsymbol e
$$
なので，運動方程式に代入してやると
$$
\boldsymbol r '' = 2 \left( \frac12 v^2 - \frac\mu r\right) \boldsymbol r - \mu \boldsymbol e
$$
ここでエネルギー
$$
E = \frac12 v^2 - \frac\mu r
$$
も当然保存量なので（束縛された楕円運動の場合はマイナスの値），最終的に
$$
\boldsymbol r'' = 2E \,\boldsymbol r - \mu \boldsymbol e
$$
というかたちが得られます．
$$
\left( \boldsymbol r - \frac{\mu}{2E} \boldsymbol e \right)'' = 2E \left( \boldsymbol r - \frac{\mu}{2E} \boldsymbol e \right) \quad (E < 0)
$$
と書くとよりはっきりと調和振動子に見えますね．
もう一回微分
天下り的に離心率ベクトルの登場していただいたのが若干悔しいので，$\boldsymbol r'''$を力ずくで計算することで運動方程式が調和振動子になっていることを再確認したいと思います．
\begin{aligned}
\boldsymbol r''' =& r \frac{d}{dt} \left( - \frac\mu r \boldsymbol r + (\boldsymbol r \cdot \boldsymbol v) \boldsymbol v \right) \\
=& -\mu \boldsymbol v + \frac{\mu \dot r}{r} \boldsymbol r + r \left( v^2 \boldsymbol v + (\boldsymbol r \cdot \ddot{\boldsymbol r}) \boldsymbol v + (\boldsymbol r \cdot \boldsymbol v) \ddot{\boldsymbol r} \right) \\
=& -\mu \boldsymbol v + \frac{\mu \dot r}{r} \boldsymbol r + rv^2 \boldsymbol v - \mu \boldsymbol v - \frac{\mu \dot r}{r} \boldsymbol r \\
=& 2r \left( \frac12 v^2 - \frac\mu r \right) \boldsymbol v \\
=& 2 E \,\boldsymbol r'
\end{aligned}
ちゃんと調和振動子のかたちになりました．
摂動も入れる
摂動も入っていないケプラー問題なんか数値計算して何の意味があるんだ（数値解法のテストにはなるけど）ということで，
$$
\ddot {\boldsymbol r} = - \frac{\mu}{r^3} \boldsymbol r + \boldsymbol f
$$
のような摂動付きのケースを考えてみます．
時間変換して，
$$
\boldsymbol r '' = -\frac{\mu}{r} \boldsymbol r + (\boldsymbol r \cdot \boldsymbol v) \boldsymbol v + r^2 \boldsymbol f
$$
もう一回微分したものには項が沢山でてきて（誰か検算してください）
$$
\boldsymbol r ''' = \bigl(2E + (\boldsymbol r \cdot \boldsymbol f)\bigr) \boldsymbol r' + 3 r^2 \dot r \,\boldsymbol f + r^3 \dot{\boldsymbol f}
$$
Hermite積分法とかを使うのならこのかたちは便利なのですが，あんまり綺麗ではありませんね．割り算が消えたので気合でもう二回ぐらいは微分できそうです