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

はじめに

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

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

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

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

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

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

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

以下で,ここの例題で使用するパッケージをロードして,警告メッセージをオフにする.
In[10]:=
Click for copyable input

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

3つの方程式からなっているが,微分項は1つしかない微分代数方程式系である.
In[11]:=
Click for copyable input

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

の導関数に指定された初期条件から開始してこの微分代数方程式系を解く.
In[12]:=
Click for copyable input
Out[12]=

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

NDSolveによって見付けられた初期条件である.
In[13]:=
Click for copyable input
Out[13]=
解のプロットを表示する.解は同じ定数値1を持つ解によって隠されてしまう.
In[15]:=
Click for copyable input
Out[15]=

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

導関数0の定常状態から始めて,解を見付けようと試みる.
In[16]:=
Click for copyable input
Out[16]=
NDSolveによって見付けられた初期条件を表示する.
In[17]:=
Click for copyable input
Out[17]=

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

方程式にを代入する.
In[18]:=
Click for copyable input
Out[18]=

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

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

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

以下はの解を求める.は微分変数なので,その値を指定することも必要である.
In[19]:=
Click for copyable input
Out[19]=
解の非零要素のプロットである.
In[21]:=
Click for copyable input
Out[21]=

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

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

NDSolveでは矛盾のない初期条件が見付けられない.
In[22]:=
Click for copyable input
Out[22]=

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

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

以下はロバートソン(Robertson)の化学反応速度問題である.大きい速度定数と小さい速度定数があるため,これはかなり硬い問題である.
In[23]:=
Click for copyable input
均衡式を微分することにより,ロバートソンの問題を常微分方程式として解く.
In[26]:=
Click for copyable input
Out[26]=

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

を表示する.
In[27]:=
Click for copyable input
Out[27]=
ロバートソンの問題を微分代数方程式として解く.
In[33]:=
Click for copyable input
Out[33]=

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

これは,常微分方程式と微分代数方程式の解の平衡式における誤差のグラフをそれぞれ黒と青で表示する.および誤差の大きさが激しく変動しているため,両対数スケールが使われる.
In[34]:=
Click for copyable input
Out[37]=

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

以下のような形式

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

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

まず振り子運動のシミュレーションの方程式,不変条件,初期条件を定義する.
In[55]:=
Click for copyable input

微分方程式は事実上不変条件の導関数であるため,その不変条件を使って方程式を解くのもひとつの方法である.

不変式を使って振り子運動の解を求める.オプションにより,NDSolveは二次方程式をについて記号的に解くのではなく,微分代数方程式として解く.
In[58]:=
Click for copyable input
Out[58]=

しかし,この解は想定されるものとは異なる可能性がある.不変式がから始まると定数という解を持つのである.実際,この初期点からでは一意の解がない.これは,について実際に解くと,関数は一意性に必要な連続性を満足しないためである.

以下では微分方程式のみを用いて振り子運動を求めてみる.ここで法が使われているのは,投影法のサブメソッドにもなるためである.
In[59]:=
Click for copyable input
Out[59]=
最後の数点上に解をプロットする.
In[60]:=
Click for copyable input
Out[60]=
NDSolveが使用した時間刻みの最後での不変量のプロットである.
In[61]:=
Click for copyable input
Out[62]=

不変量の誤差は大きくはないが,一定の割合で確実に傾斜している.これは最終的に解の信頼性に影響するほど大きくなる可能性がある.

ここでは各ステップでの運動が不変量上に存在するよう制約して,振り子運動の解を求める.
In[63]:=
Click for copyable input
Out[63]=
NDSolveが使用した時間刻みの最後の不変量を投影法でプロットする.
In[64]:=
Click for copyable input
Out[65]=
New to Mathematica? Find your learning path »
Have a question? Ask support »