前回に引き続き今回も、万有引力の法則に支配される惑星の運動を 取り上げます。問題は前回と同じです。 今回は Mathematica を使った数値計算が主題です。

微分方程式の数値解を求めるのに、NDSolve を使うのは先週説明したとおりです。 数値計算をするときは、 解くべき方程式をなるべく簡単な形にするものなので、先週 DSolveで解いた微分方程式

を NDSolve でもう一度解きましょう。アルゴリズムには Runge-Kutta 法を採用します(Runge-Kutta 法についてはここを参照)。具体的には NDSolve のオプションとして Method -> ExplicitRungeKuttaを指定します。

微分方程式を解くのは積分することと同値ですから、Runge-Kutta法 に限らず、数値積分を実行するには初期条件が必要です。DSolve で解析解を求めたとき、初期条件を指定しなくてよかった （指定しても良い。）のは、DSolve が積分定数 C[i] を自動的に 挿入して、 初期条件を指定しないことから生じる自由度を埋め合わせていた からです。しかしNDSolve は数値で解を求めなければならないので、 解を決定するには充分な初期条件を与える必要があります。 先週やった微分方程式の数値解を求めるためのアルゴリズムをきちんと理解 していれば、初期条件が決まらなければアルゴリズムのメインである逐次 計算をスタートできないことがわかるはずです。

しかし初期条件といっても、「逐次計算をはじめる時点での条件」を与える 必要はありません。それどころか初期条件として、逐次計算する範囲外での 条件を与えても構いません。Mathematica には、 「逐次計算をはじめる時点での条件」をとりあえず適当に仮定して計算を開始し、 逐次計算終了後にその解を初期条件に合うように調整する機能が実装されて います。

ここでは を初期条件として与えます。 常にθ=0での条件である必要はありません。これらの初期条件は、 解くべき方程式とともに以下のようにリスト形式で与えます。

他に DSolveと違う点は、逐次解を計算するθの範囲を指定する必要 があることぐらいです。

という出力が得られれば成功です。この出力は、 の範囲で を数値的に近似計算し、その結果を内挿（interpolation）することで、範囲 における が近似関数として求まったことを意味しています。

どんな解が求まったのか、

をお絵かきして調べることにします。それにはまず、

として、得られた解 solution を描かせたい関数に代入して、 プロット前にそれを評価しておく必要があります。これをプロットすると、

となって、θ=0 で最も重力源に近く、θ=πで最も重力源から離れており、 妥当な結果になっています。 これだけではわかりにくいので、前回と同様、カーテシアン座標系に変換して ParametricPlot で軌道を描かせてみます。 しかし NDSolve の出力、あるいはそれを代入したｒをよく見てみると リスト形式になっており、 残念ながらそのままでは ParametricPlot の入力として使うことができません。 そこで、リストからその第一要素を取り出す First を使います。実際やってみると

となって、外側にあった括弧が消えています。前にやったTake などは、リストから一部をリストとして取り出すための関数でしたから、 要素を取り出すFirst とは異なる ことに注意してください。これを使えば、

となって、前回と同じようなグラフになっていることがわかります。 例えば

とすると、

となり、前回やった厳密解と一致しています。ただし、 前回描いたのは厳密解であり、今回はあくまで数値計算による近似解 を描いていることに注意してください。つまり、 誤差が小さすぎて違いが見えなかったことになります。 ちなみに、わざと精度を落した計算をして、誤差があることを実感し てみましょう。NDSolve でオプション PrecisionGoal->3 を付けると、3桁の精度で計算することになります。

となり、4桁目から値が大きくずれていることが分かります。