Wolfram言語で解く高度な数値微分方程式入門

概要

Wolfram言語の関数NDSolveは，一般的な数値微分方程式ソルバである．この関数は常微分方程式(ODE)の大部分と偏微分方程式(PDE)の一部を扱うことができる．常微分方程式系では未知の関数 ui の数はいくつでもよいが，この関数すべてが各関数に共通の1つの「独立変数」t に依存しなければならない．偏微分方程式は2つ以上の独立変数を含む．NDSolveは微分代数方程式 (DAE)のあるものも解くことができる．微分代数方程式とは微分方程式と代数方程式が混ざったものである．

NDSolve[{eqn1,eqn2,…},u,{t,tmin,tmax}]関数 u の t についての tmin から tmax までの範囲における数値解を求める
NDSolve[{eqn1,eqn2,…},{u1,u2,…},{t,tmin,tmax}]複数の関数 ui の数値解を求める

常微分方程式の数値解を求める

NDSolveはInterpolatingFunctionオブジェクトとしての関数 ui の解を表す．InterpolatingFunctionオブジェクトは独立変数 t の値が tmin から tmax までの範囲における ui の近似を提供する．

一般に，NDSolveは反復的に解を求める．t のある特定の値から始めて，一連の段階を通り，最終的に tmin から tmax までの全範囲をカバーしようとする．

まず最初に，ui とその導関数の適切な初期条件あるいは境界条件をNDSolveに与えなければならない．これらの条件が，特定の点 t における ui[t]の値と，場合によっては導関数 ui'[t]も指定する．条件が与えられた範囲に1つの t しか存在しない場合は，等式と初期条件はひとまとめで初期値問題として扱われる．複数の点 t があると境界値問題が起る．NDSolveは記号的に標準形にすることのできる初期値問題ならほぼすべてを解くことができる（つまり，導関数の最高の次数について解くことができる）．しかし，これは線形の境界値問題に限られる．

0から2の範囲におけるtの関数xの解を，t1におけるxの初期条件を使って求める．

NDSolveを使うときの初期条件あるいは境界条件は，ui についての解が完全に得られるほど十分なものでなければならない．DSolveを使って微分方程式の記号解を求めるときに指定しなければならない制約条件の数は，これより少なくて済む．これは，DSolveが任意の記号定数C[i]を自動的に挿入することにより，明示的に指定されない初期条件に関連する自由度を表すためである．NDSolveは数値解を出さなければならないので，この種の追加的な自由度を表すことはできない．この結果，解を決定するのに必要なすべての初期条件あるいは境界条件をユーザが指定しなければならないことになる．

典型的な例として 次までの導関数を持つ微分方程式を考えてみる．この場合， 次導関数までの初期条件を与えるか，あるいは 個の点における境界条件を与えるかしなければならない．

次は二階微分方程式の初期値問題を解く．これに必要な2つの制約条件は，t0で与えられている：
次は得られた解のプロットである：
次は単純な境界値問題である：

各変数に妥当な数の条件がある限りにおいて，NDSolveを使って連立微分方程式系を解くことができる．

以下で1組の連立微分方程式の数値解を求める：
両方の解をプロットする：

初期条件はどんな種類の方程式として与えることもできる．この方程式に複数の解がある場合は，NDSolveが複数の解を生成する．

この場合の初期条件だと解が複数得られる：

平方根関数に分枝切断の問題があるために，NDSolveはy'[x]-Sqrt[y[x]^3]，y[0]-2の解を求めることができなかった．

次はNDSolveで得られた解の実数部分を示したものである（上の2つの解は厳密に実数解である）：

NDSolveは微分代数方程式(DAE)と呼ばれる，微分方程式と代数方程式が混ざった系も解くことができる．実際，上記の例はある種の微分代数方程式である．しかし，方程式は微分係数を用いて明示的に表されてはいない．典型的な微分代数方程式は微分係数について解くことは全くできず，問題は全く異なる方法で解かれなければならない．

これは単純な微分代数方程式である：

等式は両方とも微分項を持つが，変数yには微分係数がない．このため，NDSolveが警告メッセージを発している．一階の微分方程式に変換するための通常の代入を行うと，方程式のひとつが実質的に代数的になる．

また，yは代数的にしか現れないので，その値を決めるために初期条件を与える必要はない．微分代数方程式に矛盾しない初期条件を求めるのは，実際には極めて困難である．「微分代数方程式の数値解法」に関するチュートリアルでより詳細に論じてある．

これは解のプロットである：

このプロットからyの微分係数が恣意的な速さで変化する傾向があることが分かる．方程式には明示的に現れないが，この条件はこれ以上ソルバが継続できないことを意味している．

微分方程式の未知の関数は，必ずしも単一のシンボルで表さなくてもよい．未知の関数が多数ある場合は，関数にx[i]あるいはxiのような名前を与えると便利なことが多い．

25個の連立微分方程式と初期条件を構築し，これを解く：

これは実際に片端の温度が一定のロッドについての熱方程式の近似解を計算する（より正確な解が求めたければ を大きくする）．

結果は，片端で一定温度が保たれている，長さ1の一次元のロッドについての熱方程式の近似解になる．右は空間的な値とみなされる解を時間関数として示す：

未知の関数がベクトル（あるいは行列）値を持つよう指定することもできる．未知の関数の次元数はその初期条件から取られる．方程式がWolfram言語の算法規則に則った一定の次元を持つ限り，未知のスカラー方程式とベクトル方程式を混ぜても構わない．InterpolatingFunctionの結果は未知の関数と同じ次元数の値を返す．記号的に表現するのが困難である，あるいは効果的ではないプロセスが微分方程式系を支配している場合には，非スカラー変数を使うと便利である．

以下でベクトル値の未知の関数を使って上記と同じ系を解く：

より多くの独立変数を指定することで，NDSolveを用いて直接解くことができる偏微分方程式もある．

NDSolve[{eqn1,eqn2,…},u,{t,tmin,tmax},{x,xmin,xmax},…]
関数 u[t,x …]の偏微分方程式系を t が tmin から tmax まで，x が xmin から xmax, …までという風に解く
NDSolve[{eqn1,eqn2,…},{u1,u2,…},{t,tmin,tmax},{x,xmin,xmax},…]
複数の関数 ui についての偏微分方程式系を解く

偏微分方程式の数値解を求める

これはNDSolveで直接求めた熱方程式の解である．
解のプロットである：

現在NDSolveは数値的線分割法を使って偏微分方程式を解いている．このメソッドは最低でも1つの独立変数を含む初期条件を持つ問題に限られる．例えばラプラス(Laplace)方程式のような明示的な偏微分方程式は境界条件を必要とするので，この方法では解くことができない．これが解ける問題に関しては，線分割法は極めて一般的であり，偏微分方程式系や非線形系も，しばしば極めて速く解くことができる．この方法については「線の方法」で詳述する．

これは，周期的境界条件を持つ2つの空間次元まで一般化された非線形サイン・ゴードン方程式の数値解を求める：
これはt3における結果のプロットである：

上で説明したように，NDSolveは独立変数 t の連続するステップを取ることで働く．NDSolveは適応的な手続きでステップの刻み幅を決定する．一般に，解が特定の範囲で急速に変化しているような場合には，NDSolveは刻み幅を小さくして解をより正確に追跡しようとする．

NDSolveを使うと，求める解の精度や確度が指定できる．一般に，NDSolveは求められた解が指定されたAccuracyGoalかPrecisionGoalのどちらかを満たすまで徐々に刻み幅を小さくしていく．AccuracyGoalの設定値が解で許される絶対誤差を事実上決定するのに対し，PrecisionGoalは相対誤差を決定する．一般に，値が零に近付いていくような解を追跡したい場合は，AccuracyGoalの設定値を大きくする必要がある．AccuracyGoal->InfinityとするとNDSolveはPrecisionGoalだけを使う．一般に，AccuracyGoalとPrecisionGoalは特定の時間ステップの局所的な誤差を管理するために使われる．ある種の微分方程式では，この誤差が累積されることがあるため，時間間隔の終端における結果の精度あるいは確度がAccuracyGoalあるいはPrecisionGoalの設定値よりはるかに小さくなることがある．

NDSolveはWorkingPrecisionで設定した値を使って内部計算の精度を決定する．AccuracyGoalあるいはPrecisionGoalにより大きな値を設定した場合は，WorkingPrecisionの値も大きくする必要があることがある．デフォルト設定のAutomaticの場合，AccuracyGoalとPrecisionGoalは両方ともWorkingPrecisionの設定値の半分になる．

NDSolveは指定の許容値を満たしているかどうかを判定するのに誤差推定を用いる．方程式系で使うときにはオプションNormFunction->f の設定を用いて異なる要素の誤差を結合する．このノルムは許容値によってスケールされ，NDSolveが

というステップを取るようにする．ここで は誤差の 番目の要素であり， は現行の解の 番目の要素である．

微分方程式の高精度の解を生成する：
端点における解の値である：

NDSolveはその適応的な手法により，t とともにさまざまな速度で変化する複数の要素を持つ「硬い」微分方程式を解くことができる．

NDSolveは解を正確に突き止めるまで刻み幅を小さくするという一般的な手法を用いる．しかし，真の解に特異性がある場合には問題がある．この場合，NDSolveは刻み幅を永遠に小さくし続け，決して終らない可能性がある．この問題を避けるためには，オプションMaxStepsを使って解が求まるまでにNDSolveが取るステップ数の最大値を指定する．常微分方程式の場合，デフォルト設定はMaxSteps->10000である．

NDSolveは10,000ステップ後に停止する．
実際，この解にはx=0のところに特異点がある：

滑らかな解を持つほとんどの方程式では，デフォルト設定のMaxSteps->10000で十分である．しかし，解が複雑な構造の場合は，MaxStepsの設定値を大きくしなければならないことがある．設定をMaxSteps->Infinityとすると，使用するステップ数の上限がなくなる．

NDSolveには解を計算するための複数のメソッドが組み込まれており，追加的なメソッドを加えるメカニズムも備えられている．デフォルト設定のMethod->Automaticの場合，NDSolveは微分方程式に適したメソッドを選択する．例えば，方程式が硬い場合，必要に応じて陰的メソッドが使われる．方程式が微分代数方程式(DAE)であれば，DAEメソッドが使われる．一般に，微分方程式を実際に解かずにその解の性質を判定することは不可能である．従って，広範囲の問題を解くのにはデフォルトのAutomaticメソッドがよい．しかし，特定の問題に対して選ばれたメソッドが使い得る最高のものではないかもしれない．また，シンプレクティック積分器のように解の特定の性質を保存する方法を選んだ方がよいこともある．

特定の系のための適切なメソッドを選ぶことは極めて困難である．それぞれのメソッドには独自の設定があり，これが解の有効性と正確さに大きく影響するために，さらに問題は複雑になる．どのメソッドをどのような場合に使うべきか，また特定の問題に適用する際にどのようにそのメソッドを調整すればよいかの大枠が掴めるように，このドキュメントの大部分は各メソッドの説明に当ててある．さらに，NDSolveにはユーザ独自の定義を持つメソッドによる式や結果があたかも組込み関数によるもののようにNDSolveによって処理できるメカニズムが備わっている．

NDSolveは一般に3段階で解を計算する．第1段階では方程式が処理され，通常は関数の形で方程式の右辺に標準形で置かれる．第2段階では初期条件から解を反復させるのに，この関数が使われる．最終段階では反復処理中に保存されたデータが処理されてひとつあるいは複数のInterpolatingFunctionオブジェクトになる．関数をNDSolve`コンテキストで使うとこれらのステップを別々に行うことができ，さらに重要なことに，反復処理がコントロールしやすくなる．各ステップは微分方程式を解くのに必要なすべてのデータを保持しているNDSolve`StateDataオブジェクトによって繋がれている．