惑星の運動 - 2

DSolve は賢い関数で、リストであらわされる連立はもとより、線形、非線形を問わず常微分方程式の一般解を見つけることができます。しかしこれまで見てきたように、解ける場合のほうが実は限られていて、あまり実用的ではありません。

これを改善することを考えましょう。ほとんどの人は解析力学を習っているでしょうから、そこで Newton の運動方程式に変わって登場する正準方程式を知っているはずです。これは1階の微分方程式ですから、対応する Newton 方程式より数学的な取り扱いが楽なはずです。

c)正準方程式を解く

同じ問題を正準形式に直すと、正準方程式は以下のようになります。 (解くべき微分方程式の導出方法については、ここを参照してください。[ 訂正：２つ目の式で、V= -(1/mu)(1/r) --> V= - mu (1/r)])

q₁=r、q₂=θ と一般化座標を導入し p₁、p₂をそれぞれの一般化運動量としました。

これを Mathematica に解かせるには、

とします。以上で解けるはずなのですが、試してみるとまた意味不明の表式が出てきます（本当は意味があるのですが）。やはりこの3つの微分方程式を手動で連立させてから解かせないと無理なようです。ちなみに手動で連立させると、簡単に先ほどの最終的な微分方程式に帰着させることができるので各自やってみてください。

このように正準方程式を用いると、先ほどよりはかなり見通しがよくなるものの、やはり限界が見えてきてしまいます。そこで登場するのが数値的に微分方程式を解こうという方法です。

d)微分方程式の数値解（アルゴリズム）

一般に、解析解は限られた（特殊な）条件下でのみ厳密解が得られますが、わずかに拡張しようとするととたんに厳密解を求めることは困難になります。このような場合、Mathematica を使って数値解を求めることで、解ける対象を広げることができます。詳しくは次回の演習で行いますが、今回の演習ではその準備として、計算のアルゴリズムや「NDSolve」というコマンドの使い方を紹介します。

数値計算は奥が深く、完全に計算機に任せていると危険なこともあります。そのため、数値計算の原理を理解しておくことが必要不可欠です。そのためまずここでは、微分方程式の数値解を求めるアルゴリズムを解説します。といっても、このトピックだけで本が１冊かけるほどなので、ほんの入門程度の内容しか扱いません。興味がある人は数値計算に関する教科書を参照してください。

アルゴリズムへ (一通り読み終えたら「BACK」をクリックして戻ってきてください)

e)NDSolveの練習

さて、アルゴリズムの概要が理解出来たら、少し実践してみましょう。以下ではNDSolveという微分方程式を数値的に解くコマンドを扱います。今回の演習では準備として、このコマンドをいくつかの方程式に適用してみることにしましょう。

s = NDSolve[{y'[x] == Sin[x], y[0] == 1}, y, {x, 0, 30}]

これを実行するとy→ Interpolating...のようなよくわからない表示が出ると思いますが、これは数値的にちゃんと解が得られていることを示しています。次に得られた解を以下のコマンドによりプロットしてみましょう。

Plot[Evaluate[y[x] /. s], {x, 0, 30}]

このコマンドの意味は次の通りです。最初の計算で、sの値には「y→得られた解」という「置き換えのルール」が代入されてます。（Mathematicaではこのように、変数の中に「ルール」を格納しておくこともできるのです。） Evaluate[y[x] /. s]の部分は、yにこのルールを適用してから、その値を評価せよという意味です。そして（少々回りくどいですが）その評価したものを最終的にプロットしているわけです。今回の微分方程式は手でも解けるものですから、自分の手で解いてみて、それがプロットの結果とあっているか確認してみましょう。また、以下の例題にもチャレンジしてみましょう。

例題１:手では解けない場合

例題２：2変数の場合

他にも自分で解けるor解けない方程式を書いてみてプロットしてみましょう。また、境界条件を満たす解が存在しない場合どのように表示されるか確かめてみましょう。