MATHEMATICAチュートリアル

微分代数方程式の数値解法

はじめに

一般に，常微分方程式系(ODE)は以下のような標準形で表される．

従属変数の導関数は，独立変数と従属変数で明示的に表される．関数が十分な連続性を持つ限り，独立変数の特定の値について従属変数の値が与えられる初期値問題で，常に一意解が求まる．

微分代数方程式(DAE)では，導関数は一般に明示的に表されない．事実，従属変数の導関数には，通常方程式の中に含まれないものもある．微分代数方程式系の一般形は

でる．ここで，のヤコビアンは特異である可能性がある．

微分代数方程式系は，独立変数で微分することにより常微分方程式系に変換できる．微分代数方程式の「階数」は，実質的には常微分方程式系を得るために微分代数方程式を微分する回数である．たとえ微分が可能であっても，微分方程式の数値シミュレーションによってもとの微分代数方程式の特性が失われることがよくあるので，微分は計算方法としては通常使用されない．

従って微分代数方程式の数値的メソッドは微分代数方程式系の標準形で使用できるよう設計されている．NDSolveのメソッドは通常一階微分代数方程式を解くよう設計されているが，これより大きい階数の問題でも使えることがある．

このチュートリアルでは，微分代数方程式系と常微分方程式系の解法の違いのいくつかを示す多数の例題を提示する．

以下で，ここの例題で使用するパッケージをロードして，警告メッセージをオフにする．

In[10]:=

初期条件の指定は常微分方程式と微分代数方程式とでは全く異なっている．前述の通り常微分方程式では，初期条件の集合によって解が一意に決定される．微分代数方程式ではこれほど単純にはいかない．方程式を満足する初期条件を見付けることですら困難な場合があり得る．この問題をよく理解するために，以下に例を示す[AP98]．

3つの方程式からなっているが，微分項は1つしかない微分代数方程式系である．

In[11]:=

初期条件は明らかに自由ではない．2番目の方程式でが0か1になることを要求している．

の導関数に指定された初期条件から開始してこの微分代数方程式系を解く．

In[12]:=

Out[12]=

この解を得るために，NDSolveはまず関数SolveにFindRootのような手続きを組み合せて，方程式を満足する初期条件を検索する．矛盾のない初期条件が見付かったら微分代数方程式法で微分代数方程式が解かれる．

NDSolveによって見付けられた初期条件である．

In[13]:=

Out[13]=

解のプロットを表示する．解

は同じ定数値1を持つ解

によって隠されてしまう．

In[15]:=

Out[15]=

しかし，方程式を満足する初期条件すべての解があるとは限らない可能性がある．

導関数0の定常状態から始めて，解

を見付けようと試みる．

In[16]:=

Out[16]=

NDSolveによって見付けられた初期条件を表示する．

In[17]:=

Out[17]=

が1に設定された方程式を見ると，以上に進めない理由が分かる．

方程式に

を代入する．

In[18]:=

Out[18]=

真ん中の方程式が実質的になくなった．最後の方程式をで微分すると，条件が得られるが，最初の方程式が初期条件の値と矛盾してしまう．

の唯一の解はであり，この解によってこの系の階数は2であることが分かる．

この問題の別の解集合はのときにある．これを初期条件に指定することで解を求めることができる．

以下は

の解を求める．

は微分変数なので，その値を指定することも必要である．

In[19]:=

Out[19]=

解の非零要素のプロットである．

In[21]:=

Out[21]=

通常，パラメータ化された一般解が存在するので，微分変数の初期条件を指定しなければならない．このの問題では，一般解はなので，解を決定するためにはを与える必要がある．

方程式に矛盾しない初期条件を見付けることは難しい問題となることがあるので，NDSolveで常にそれを見付けることができるとは限らない．「多くの場合，アプリケーションで微分代数方程式系を解くときに最も困難な部分は，計算開始に使用する矛盾のない初期条件の集合を決定することである．」[BCP89]

NDSolveでは矛盾のない初期条件が見付けられない．

In[22]:=

Out[22]=

NDSolveが矛盾のない初期条件を見付けることができない場合，よい初期値でFindRootを使うか，別の手続きを使って矛盾のない初期条件を得て，それを使うかすることができる．よい初期値に近い値が分かっているなら，NDSolveはその値を使って検索を始めるので，役に立つ場合がある．従属変数とその導関数の値を指定することもできる．

一階の微分代数方程式系では，常微分方程式ソルバを微分・利用して解が得られることがよくある．

以下はロバートソン(Robertson)の化学反応速度問題である．大きい速度定数と小さい速度定数があるため，これはかなり硬い問題である．

In[23]:=

均衡式を微分することにより，ロバートソンの問題を常微分方程式として解く．

In[26]:=

Out[26]=

問題の系が硬いことは，との主要な変動が2つの完全に異なる時間スケール上にあることにより裏付けられる．

解

と

を表示する．

In[27]:=

Out[27]=

ロバートソンの問題を微分代数方程式として解く．

In[33]:=

Out[33]=

与えられた要素についての解は非常に近く見えるが，化学平衡の条件を比較すると，その違いが分かる．

これは，常微分方程式と微分代数方程式の解の平衡式における誤差のグラフをそれぞれ黒と青で表示する．

および誤差の大きさが激しく変動しているため，両対数スケールが使われる．

In[34]:=

Out[37]=

この場合，解は両方とも想定される許容誤差を大きく上回って平衡式を満足した．微分代数方程式の解の平衡式における誤差が大きくなっている点があっても，全体を見た場合，突然の変動区間が過ぎると微分代数方程式の解は制約条件を満足するよう引き戻されているという点に注目されたい．

以下のような形式

の微分代数方程式を解きたい場合があるとしよう．ここで，微分方程式の解は特殊な制約条件を満足しなければならない．この微分代数方程式は階数が大きすぎるし，NDSolveでは方程式の数が従属変数の数と同じであると想定されるので，このような微分代数方程式をNDSolveで直接扱うことはできない．その上，しかし，NDSolveにはProjection法があり，これを使ってこの問題を解く場合がある．

このような制約条件付きの系の非常に簡単な例として，振り子運動をモデル化した非線形振動がある．

まず振り子運動のシミュレーションの方程式，不変条件，初期条件を定義する．

In[55]:=