ケプラー問題の数値解


逆二乗の中心力のもとで物体がどのように運動するかを調べるケプラー問題。微分方程式を解く能力がなければ楽しめないかと言えば、そういうわけでもない(物理を勉強する人にはもちろん必須)。

最近の強力なプログラミング環境を使えば、物体の軌道を簡単に視覚化できる。インタラクティブなものにするのも簡単(さすがに、微分方程式は立てられないといけない)。

Mathematicaならこんな感じ(インタラクティブに操作できるのは、重力の強さを決めるkと時間の上限tmax、初速度v0、初期座標pos)。Wolfram CDF Playerがインストールされていれば実際に動かすことができる(UMMでも動く)。

Manipulate[
 s = NDSolve[{
    x''[t] == -k x[t]/(x[t]^2 + y[t]^2)^(3/2),
    y''[t] == -k y[t]/(x[t]^2 + y[t]^2)^(3/2),
    x[0] == pos[[1]], y[0] == pos[[2]],
    x'[0] == v0[[1]], y'[0] == v0[[2]]},
   {x[t], y[t]},
   {t, 0, tmax}];
 ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, tmax}],
 
 {{k, 5}, 1, 10},
 {{tmax, 10}, 1, 100},
 {{pos, {2, 1}}, Locator},
 {{v0, {-1, -1}}, {-5, -5}, {5, 5}}]


Mathematicaの微分方程式ソルバーNDSolveはかなり強力で、このエントリのような単純なものだけではなく、摂動があるような問題も扱える。

解析的に解きたい人は、ランダウ=リフシッツ『力学・場の理論』のような教科書を参照して欲しい(福島登志夫『天体の位置と運動』の方がやさしいかな )。

参考:Mathematicaで解くケプラー問題

ジンマーマン, オルネス『物理学のためのMathematica』では、ハミルトン=ヤコビ方程式を使う方法が紹介されている。

関連:プリンキピアを読む

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です