|
前回に引き続き今回も、万有引力の法則に支配される惑星の運動を 取り上げます。問題は前回と同じです。 今回は Mathematica を使った数値計算が主題です。 |
| a) 微分方程式の数値解 -- NDSolve |
|
微分方程式の数値解を求めるのに、NDSolve を使うのは先週説明 したとおりです。数値計算をするときは解くべき方程式をなるべくシンプル な形にするものなので、先週 DSolveで解いた微分方程式
を NDSolveでもう一度解きましょう。 アルゴリズムには Explicit Runge-Kutta法を採用します。具体的には NDSolve のオプションとして Method -> "ExplicitRungeKutta" を指定します。 微分方程式を解くのは積分することと同値ですから、Explicit Runge-Kutta法 に限らず、数値積分を実行するには初期条件が必要です。DSolve で解析解を求めたとき、初期条件を指定しなくてよかった （指定しても良い。）のは、DSolve が積分定数 C[i] を自動的に 挿入して、 初期条件を指定しないことから生じる自由度を埋め合わせていた からです。しかしNDSolve は数値で解を求めなければならないので、 解を決定するには充分な初期条件を与える必要があります。 先週やった微分方程式の数値解を求めるためのアルゴリズムをきちんと理解 していれば、初期条件が決まらなければアルゴリズムのメインである逐次 計算をスタートできないことがわかるはずです。 しかし初期条件といっても、「逐次計算をはじめる時点での条件」を与える 必要はありません。それどころか初期条件として、逐次計算する範囲外での 条件を与えても構いません。Mathematica には、 「逐次計算をはじめる時点での条件」をとりあえず適当に仮定して計算を開始し、 逐次計算終了後にその解を初期条件に合うように調整する機能が実装されて います。
ここでは
[訂正：Method -> ExplicitRungeKutta を Method -> "ExplicitRungeKutta" として下さい。]
という出力が得られれば成功です。この出力は、
|
| b)惑星軌道の視覚化 |
|
どんな解が求まったのか、
をお絵かきして調べることにします。それにはまず、
として、得られた解 solution を描かせたい関数に代入して、 プロット前にそれを評価しておく必要があります。これをプロットすると、
となって、θ=0 で最も重力源に近く、θ=πで最も重力源から離れており、 妥当な結果になっています。 これだけではわかりにくいので、前回と同様、カーテシアン座標系に変換して ParametricPlot で軌道を描かせてみます。 しかし NDSolve の出力、あるいはそれを代入したｒをよく見てみると リスト形式になっており、 残念ながらそのままでは ParametricPlot の入力として使うことができません。 そこで、リストからその第一要素を取り出す First を使います。実際やってみると
となって、外側にあった括弧が消えています。前にやったTake などは、リストから一部をリストとして取り出すための関数でしたから、 要素を取り出すFirst とは異なる ことに注意してください。これを使えば、
となって、前回と同じようなグラフになっていることがわかります。 ただし、前回描いたのは厳密解であり、今回はあくまで数値計算に よる近似解を描いていることに注意してください。例えば
とすると、前回やった厳密解であれば
となったのですが、今回はどうなっているでしょう？ ちなみに、NDSolve でオプション Method -> "ExplicitRungeKutta" を指定しないで計算するとどうなるでしょうか？ 省略時にはアルゴリズムとして多重ステップ法が用いられますから、 得られる計算結果にも微妙な違いが生じるかもしれません。 各自、確かめてみてください。 その他、オプションを指定することで、様々な方法で微分方程式を（近似的に）解くことが できます。時間があれば試してみて下さい。 |
を初期条件として与えます。
常にθ=0での条件である必要はありません。これらの初期条件は、
解くべき方程式とともに以下のようにリスト形式で与えます。
の範囲で
を数値的に近似計算し、その結果を内挿（interpolation）することで、範囲
