2014年11月9日日曜日

Rで微分方程式

微分方程式を解かず(積分せず)に&楽に解(結果)を知りたい。大学で習うのは簡単に積分できる簡単な微分方程式ばかりですが、興味深い微分方程式の多くはそうした簡単に積分できるものばかりではありません。むしろ興味深い性質を持つのは簡単に積分できないからということも多く、そうした場合にどうするかがいつも頭の痛い問題でした。もちろんMathematicaを用いるとか、お金がないならOctaveを使うとか、微分方程式をコンピュータに解かせることはできるのでしょうが、それを学生にやらせるにはいろいろな問題があります。例えば、フリーでないので、学生の自習用に同じ環境を学生のパソコンに用意するのが難しいとか、数値計算でしか使わないソフトを習得させるためのコストを学生は払う価値があるのか、なのです。

そこで、フリーで自由にインストールして自習環境が容易に準備可能&数値計算専用ソフトと遜色ない数値解析能力&他の用途にも使用可能という3つの条件を満たすものとして、RとdeSolveというパッケージに着目し、3年前期の選択科目「コンピュータ数値計法」で扱うことにしました。2年後半の「統計データ解析」でRを使ってデータの統計解析(記述統計、回帰分析、母数の推定・検定など)を扱っているので、同じRの環境で数値計算を学んだほうが学生のコストを抑えられるという考えもあります。

deSolveは微分方程式を数値積分するRのパッケージで2010年に公開されたもの。もともとRにはodesolveという常微分方程式を数値積分するパッケージ(2001)があったのですが、それを徹底的に強化し、常微分方程式、拘束条件のある常微分方程式から、遅延のある常微分方程式、偏微分方程式を数値積分する環境としたものです。


このdeSolveのチュートリアルとして書かれたのが、K.Soetaert,J.Cash and F. MazziaによるSolving Differential Equations in Rというテキスト。K.SoetaertはdeSolveを開発者の一人です。



テキストは数値積分のアルゴリズムとdeSolveで用いるアルゴリズム、計算誤差の議論から始まり、さまざまな微分方程式をdeSolveで解く方法が解説されています。

例えば、左に示した二つの図は、第8章の偏微分方程式の差分解法のあとの第9章で扱っている熱伝導方程式、波動方程式、ラプラス方程式、1次元および2次元のブッルッセレータ方程式、サインゴルドン方程式、非線形シュレーディンガー方程式のうちの2次元ブルッセレーター方程式(上の図)と非線形シュレーディンガー方程式(下の図)を数値積分した結果を示したものです。

ソースコードは初期条件や格子の差分化の指定などを含めても10数行の短いもの。それでこうした結果を得られるのは素晴らしいの一言。時間が出来たら2次元ナウ゛ィエストークス方程式を数値積分して対流から乱流になる様子をローレンツ方程式と比較してみたいと思います。示したいのは波長の短い撹乱が長波長側には影響せず、バタフライ効果がないこと。つまり、ブラジルで蝶が羽ばたいてニューヨークでハリケーンが起きるなんてバカなことはないことを数値計算で検証したいだけなのですが、deSolveを使えば楽に出来そうで楽しみです。

下は2次元サインゴルドン方程式の解のスナップショット。