Excelで解くケプラーの方程式

「黒点の移動速度の求め方」で太陽の自転軸の傾いている方向を表す式として

　　θi = θ + π(ti - t1)/1[year]

と書きました。

θは太陽の「自転軸の傾いている方向」なのですが、これは太陽から見た地球の方位角（の逆方向）（に一定値を足したもの）でもあります。

上の記事は「地球の軌道は円軌道」という前提ですからこれでいいのですが、じっさいの地球の軌道は楕円なのでこれは正確ではありません。

正しくθに相当する平均近点角lに対する太陽から見た地球の方位角f（真近点角）を求める必要があります。

それにはケプラーの方程式

　　u - e sin u = l

を解いて平均近点角lから離心近点角 u を求める必要があります（e:離心率です。真近点角fは離心近点角uから求めることができます）

この方程式は（非線形で）解析的には解けないので数値的に解を求める必要があります。

離心近点角uの初期値を適当に決めて上式の残差

　　R= l - u + e sin u

を求め、残差Rが0に（近く）なるように計算を進めていきます。

残差 R を u で微分すると

　　dR/du = -1 + e cos u

となりますから、まず

　　R1 = l - u1 + e sin u1

を求め

　　Δu1 = - R1 / ( dR/du ) = - R1 / ( -1 + e cos u1 )
　　．
　　．
　　Δui = - Ri / ( dR/du ) = - Ri / ( -1 + e cos ui )

から

　　u2 = u1 + Δu1
　　．
　　．
　　ui = ui-1 + Δui-1

として新たな R を求めるということを繰り返しRが十分小さくなるu2,..,uiを求めればいいはずです。

u1をどう決めるか問題ですが

　　平均近点角≒離心近点角

と考え

　　u1 = l + e sin l

とします。

Excelでこれを計算してみます。

各セルに次のように計算式をセットします。
離心率eは0.01672とします。これは地球のじっさいの離心率です。
平均近点角l=1（57.29...度）の場合を計算します。


じっさいに計算するとあっけなく収束してしまい離心近点角 u=1.01420 が求まります。


-----

ところでExcelにはソルバーという機能があります。これを使うともっと安直です。

「ソルバー」はアドインなので最初にインストールする必要があるのですが、一度インストールしてしまえばExcelの一つの機能としてふつうに使えます。

次のようなシートを作ります。


データ=>分析=>ソルバーからソルバーを起動します。

uつまりA5を変化させてRつまりB5セルを0にすることが目標ですから
　　目的セル： ＄B$5
　　目標値： 値： 0
　　変化させるセル： ＄A$5
と設定します。


そして実行ボタンをクリックすると次の結果が一瞬で得られます。


　　u=1.01420

となりもちろん上の繰り返し計算を行ったときと同じ結果が得られます。

------

「Excelで解くケプラーの方程式」　編集