常微分方程式にトライ（差分法） - 理系ツールのクールな使いこなし

BB非解析解

1.　数値解析解

AAでは解が陽的に数式で得られるケースを挙げたが、自然界の振舞を表すODEはそれだけではない。非線形現象のためむしろ数式で表せない方が圧倒的に多い。そんな時は数値解析手法に頼る事になる。つまりAAで頑張っても解けない場合は数値解析しましょうという訳である。このブログの性格上、コーディングは最終手段なので、できる限りソフトウェアの搭載機能を利用していきたい。

数値解析手法は大きく次の３つに分けられる。陽解法は数値解析の基礎である。とにかくODEをどのようにコンピュータに計算させるのかを理解するために通らなければならない入口だ。本ブログではルンゲクッタ法を使った計算例に限定する。陰解法は非線型ODEがターゲットである。ODEの性質により著しく解きにくい（stiff）場合などは予測子—修正子法が使用される。

（i）陽解法

・オイラー法：微分方程式の数値解法の一番手は、どのテキストを眺めてもオイラー法から始まる。解りやすいのが一番だ。しかし一次精度のため計算刻みを相当ゼロに近づけないと計算誤差が累積していき、近似解精度が保てなくなる。理解しやすいがコスパが悪すぎる。

・ルンゲクッタ法：2次は修正オイラー法とも言われる（2次のテイラー展開利用）。一般的なのは4次のテイラー展開を利用したもので、よくRK4と表示される。5次のテイラー展開利用したものも存在する。当然次数の高い手法ほど時間がかかるが精度は高くなる。

　(ii) 陰解法

・後退オイラー法：一次精度なので、通常は使用する必要はない。

・クランク・ニコルソン法：二次精度。

　(iii) 予測子—修正子法：予測子に陽解法を、修正子に陰解法を用いる。計算回数が多くなると振動を起こす。

・Milne

・Adams-Moulton法

計算例

　（1）数式処理ツールMapleの計算結果を以下に示す。

大きく二つに分けてある。一つはグラフ的解法。withコマンドでパッケージDetoolsをロードして、DEplotを走らせると解の流れが見られる。計算しなくてもODEの解の全体像が俯瞰できる、凄いよね？　初期条件を入力することで、当然それを通過した解曲線がプロットされる。ここで大事なことは、初期条件により解の結果が異なる、つまり初期値に敏感なODEが存在する事を認識することだ。

二つ目はMapleに搭載された様々なルンゲクッタ法を使用した例である。rk3、rk4、rkf45、ck45、dverk78、rosenbrock。計算結果はリスト表示しているが、その違いは判別できないくらい小さい。コーディングする前に色々試せるということ。

------------------------------------------------------------------------------------

図1

------------------------------------------------------------------------------------

（2）次はWolfram CloudのMathematicaの計算結果だ。NDSolveコマンドでODEの数値解析ができる。計算結果はオブジェクトInterpolatingFunctionに格納される。そのデータをEvaluateコマンドを使用してplotする。

二つ目はMethodにExplicitRungeKuttaを指定したケースである。マニュアルにはImplicitRungeKuttaも記載されている。

両プロットとも違いは分からないほどだ。

------------------------------------------------------------------------------------

図2

--------------------------------------------------------------------------------------

（3）最後はLiveMathによる計算結果だ。

数値解析する場合はまず計算したいODE全体をセレクトする。画面上部の計算メニューからTabular->微分方程式の数値解法をセレクトする。すると図4の様なODEソルバー画面が現れるので、画面上部のボックスからルンゲ-クッタ法をセレクトする。また独立変数の範囲や計算間隔を入力すると、計算結果の個数を合計欄に表示してくる。間隔を小さくすれば、それだけ計算合計が増えるのでプロットは綺麗になる（時間は掛かる）。初期値は初期値画面で入力するが、従属変数の値しか入れられない（独立変数の初期値は範囲の最初の値になる）。okボタンを押すと計算結果が表（図3のプロットのすぐ上にある行）で格納される。そのデータをプロットしたものが図3である。

-------------------------------------------------------------------------------------

図3

図4

図5

-------------------------------------------------------------------------------------

今回の3種類のソフトウェアはコマンドこそ異なるが、計算結果が非常に近いことが理解できるはずだ。

実は、AAの解析解グループにはグリーン関数法・たたみ込み積分法などの手法も存在するが、どのような場合にこれらを使用しなければならないのか、不勉強で残念ながら僕では伝える力が無い。また数値解析の前段階には近似解法として、テイラー級数やフーリエ級数を利用した手法もある。詳細は参考書を挙げるだけで勘弁願いたい。

グリーン関数法・たたみ込み積分法は