今回は、万有引力の法則に支配される惑星の運動を取り上げます。 2 体の運動を考え、一つを太陽、一つを惑星と見なします。 万有引力は中心力、 すなわち角度方向には力が働かないという特徴がありますから、 太陽を中心とし、惑星の軌道面を含む 2 次元極座標で考えるのが簡単です。

Mathematica を使って常微分方程式を解く方法は前にも取り上げましたが、 今回はその復習とその限界の認識、そしてそれを乗り越えるための方法 （微分方程式の数値解法）の基礎を勉強します。

Newton の運動方程式は

です。ベクトルrは太陽から惑星への位置ベクトルです。 μ = GMm です (G: Newton の万有引力定数, M :太陽の質量, m :惑星の質量)。 2 次元極座標を使うと下の連立微分方程式が得られます。

これらの微分方程式の導出についての詳細は、 ここを見て ください。

振動の視覚化では微分方程式を解くのに DSolve を使いました。 そのときもコメントしましたが、この関数はリストを使って、 今回扱うような連立微分方程式を解くこともできます。 早速試してみましょう。

eq1 は動径方向の微分方程式、eq2 は角度方向の微分方程式です。 これを解かせると

となって、計算してくれないようです。 解析解が存在しない場合には、 Mathematica はその旨の警告メッセージを出力するのですが、 どうも微分方程式が複雑すぎるのか、それすらやってくれていない ようです。しかし、 eq2 だけなら解くことができます。

記号 K$87 は積分変数を表していますが、 数字は変わる場合があります。

完全に一般的に解いてしまって、θのt依存性が求まっています。 実はC[1] は(単位質量当たりの)角運動量を表しています (詳しくはここを見てください)。 完全に解いてしまったθは用いずに、dθ/d t = (L/m)/r² (L は角運動量)を使って eq1 を解いてみましょう。

要するに、解析解が得られませんでしたというメッセージを出して 複雑な式を出力しました。しかし、 これは人間には分かりづらい形なので、 解くべき微分方程式をもう少し簡単にすることを試みた方が見通しがよさそうです。 簡単な考察より最終的な微分方程式の形が以下のように得られます。

ここで、u(θ) = 1/r(θ) です。 この方程式からは惑星のある時刻の位置が得られるのではなく、 惑星の軌道が得られます(rとθの関係式)。 これを Mathematica に解かせてみます。簡単のため

に選ぶと、

となり、やっと一般解が得られました。 任意定数は初期時刻での位置と速度を与えれば求まります。たとえば、 惑星の初期位置の角度は座標軸の取り替えでいつでもゼロにできるので、 そのようにとると、 u(θ(t=0))=u(θ=0)=1/r(t=0) からC[1]は初期位置か決まり 一方上の式を時間で一回微分して時刻をゼロとおけば C[2] は初期速度で書けていることが分かります。

どんな解が求まったのか、お絵かきして調べてみます。 軌道の形を見るためにr(θ)をθ=0 から 2 πまで動かし、 視覚化しやすいように前回と同様、カーテシアン座標系に変換します。 三角関数の加法定理よりC[1]Cos[θ]+ C[2] Sin[θ] = A Cos[θ+ B] の形にいつでも変形できますから、θ=0 から 2 π まで動かすため B の値は軌道の形そのものには影響を与えません (Bを変えると軌道の軸の傾きが変わります)。 そこでとりあえず B = 0 とします。 残った任意定数Aの意味は視覚化のあとに説明することにして、 とりあえずA=εと置いておきます。

;

ε=0.8, 1.0, 1.2 の場合を描いてみましょう。

原点が太陽に相当します。ε < 1 の時は楕円(特にε=0で円)、 ε=1の時放物線、 1 < ε の時は双曲線になります。 双曲線の時、図に直線が斜めに横切って出ていますがそれは 実在しません(Mathematica が距離が無限大の所にある点同士を結 んでしまっているのでこのようになってしまいました)。 あと、双曲線には二つ曲線がでていますが、実在するのは太陽(原点) をくるっと回り込んでいっている方です。

このように、解はεの大きさによって、楕円、放物線、双曲線を描きます。 これらの2次曲線は円錐曲線と呼ばれ、次のような式で表されます。

εは離心率とか焦点率などと呼ばれる正の定数で、楕円の場合 0<ε<1 になります。長軸の長さは 2a になります。