科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

ケプラー方程式をニュートン法で解くときの開始点の選び方は、悩みが尽きないのです

ケプラー方程式

 M = U - e*sin(U)

を解くとき、よくニュートン法が使われます。 離心率が大きい長楕円軌道の場合、開始点の選択を誤ると繰り返し計算がバタついて、 収束しないことがあります。 ケプラー方程式は、下に凸な関数なので、 大きい側から近づくと動作が安定します(作図すれば明らかですね?)。

とにかく安定させたいなら

範囲を 0~π に限定したうえで、開始点としてπを選べば、安定した動作が期待できます。 繰り返し回数が多くなりそうですが、安定した動作には代えられません。 ニュートン法は二次収束するので、開始点が離れていても、すみやかに正解に到達します。

もう少し正解に近い初期値が欲しいなら

平均近点角と離心近点角の差に外接する3つの直線を使って、 正解より少し大きい初期値を求めることができます。

 U0 = min(\frac{M}{1-e}, M+e, \frac{M+πe}{1+e})

これは、Nijenhuis の初期値と呼ばれます。 式の意味が幾何学的に明確なところがお勧めポイントです。

近点角が小さいところの相対誤差が気になるなら

Nijenhuis の初期値が万能かというとそうでもなくて、e が 1 に近くて(例えば e=0.9 以上) M がゼロに近いところだと、初期値の相対誤差が大きくなります。 この誤差を抑えるのに、以下の式が役に立ちます。

 M + (6M)^{1/3}

離心率 e = 1 を仮定した近似解をベースにしています。実際の離心率を反映 すれば、さらに改善が見込めます。ただし、近似解は正解より大きい値を返す 必要があります(イタレーションの安定な収束を保証するため)。初期値の計算に三乗根を 使うのはやりすぎの感もありますが、全区間に渡って相対誤差を一定以下に抑えることができる点は魅力的です。

f:id:sken20k:20180313215559j:plain  U3:nijenhuis U4:三乗根

M ゼロに近い領域の挙動を確実におさえておきたいならば、 M を対数でサンプリングして、上のグラフを相対誤差でプロットすると良いと思います。

収束判定の一工夫

上で紹介した初期値を使えば、大きい側から正解に近づきます。 そうはいっても数値誤差の関係で修正方向が反転することもあります。 そういう場合に備えて、「解の振動を検出したらそこで停止する」のも一案です。

python の デコレーターを使ってみよう !

定義域(値域)を 0~π に限定した関数を実装したあと、 定義域(値域)を-π~πに拡張することを考えます。 python なら、デコレーターを使って、ちょっとおしゃれに実装できます。

def deco_expand_sign(f):
    # 符号拡張
    def _f(x):
        if x < 0:
            y = -f(-x)
        else:
            y = f(x)
    return _f

@deco_expand_sign
def keq_eli(mu, m):
    # m
    ...
    return u
    

deco_expand_sign は、双曲線にも放物線にも適用できるあたりががミソです。 楕円軌道の場合は、+/-n 周回の不確定性を解決するデコレーターの採用も見通しが良いですね。 私的には、戻り値に +/-n 周回の情報を残した設計(定義域も値域も -∞~∞)が好みです。

ケプラー方程式の解法に関しては、天文台の福島登志夫さんによる1980年代の仕事が完成度高いです。 そのぶん、実装も検証も片手間で...とはいきませんので、採用するにはそれなりの覚悟がいります。