NDSolveの"Projection"法

はじめに

微分方程式系がある種の構造を持つ場合,数値積分法がその構造を保存すると利点が多い.ある場合は,解が制約される微分方程式を解くと役に立つことがある.投影法は,数値積分法で時間刻みを取り,実際の解が発展していく多様体上に近似解を投影することで作用する.

NDSolveには,適切であろうと思われる微分代数ソルバが含まれている.これについては「微分代数方程式の数値解法」でより詳細に記述してある.

方程式の形式は,DAEソルバによって要求される形式に還元できないことがある.さらに,いわゆる指標低減法では,微分方程式系が持ち得るシンプレクティック性等の構造特性を破壊する可能性がある([HW96]と[HLW02]を参照のこと).この例は,DAEのドキュメントに記載してある.

このような場合,制約条件が確実に保存できるように,微分方程式系を解いてから投影法を使うことができる場合が多い.これが"Projection"メソッドの考え方である.

微分方程式系が について可逆であれば,対称投影法が有利である([H00]を参照のこと).対称投影法は通常の投影よりも一般にコストがかかるため,NDSolveにはまだ実装されていない.

不変量

以下の微分方程式

で, はべクトルあるいは行列である.

定義:非定数関数 は,すべての について が成り立つとき,(1)の不変量と呼ばれる.

これは,(2)の解 すべてがを満足するということである.

不変量と同義語の「第1積分」,「保存量」,「運動定数」等もよく使われる.

多様体

である次元部分多様体があるとする.

微分方程式(3)があるなら, はすべての について ということである.これは不変性よりも弱い仮定であり,は弱い不変量といわれる([HLW02]を参照のこと).

投影アルゴリズム

が1ステップの数値積分の解を表すとする.制約条件付きの最小化問題は,以下の系を導く([AP91],[HW96],[HLW02]を参照).

作業を保存するために,と近似される.(4)の最初の式を2つ目の式に代入すると, について以下のように簡略化されたニュートン法が導かれる.

のとき,

最初の増分の大きさはなので,(5)は通常迅速に収束する.

高次の積分法を使うときに追加されるコストは,投影ステップでニュートン法の反復回数を少なくすることで相殺できる.

"Projection"メソッドの終了基準には,オプション"IterationSafetyFactor"NDSolveによって使われる作業精度の1ULP(Unit in the Last Place)と結合される.

例題

ユーティリティパッケージをいくつかロードする:

一次不変量

化学反応をモデル化する硬い系を定義する:
この系は線形不変量を持つ:
一次不変量は一般に,その不変量の誤差のプロットで観察されるように,デフォルトのNDSolveメソッドをはじめとする数値積分([S86]を参照)によって保存される:

従って,この例では"Projection"メソッドを使う必要はない.

数値法の中には,正確に二次不変量を保存するものもある([C87]の例を参照のこと).このようなメソッドのひとつに,陰的中点公式(あるいは1段のガウス公式,陰的ルンゲ・クッタ法)がある.

調和振動

調和振動を定義する:
この調和振動は以下の不変量を持つ:
この系を"ExplicitRungeKutta"メソッドを使って解く.不変量の誤差が,ほぼ線形に増加する.これはハミルトン系に適用された散逸的なメソッドでは通常の動作である:
これも"ExplicitRungeKutta"メソッドを使って系を解くが,これは各ステップの最後に解を投影する.不変量の誤差のプロットにより,解が四捨五入するまで保存されていることが分かる:
これはハミルトン系(不変量がハミルトニアン)なので,シンプテクティック積分はこの問題にうまく作用し,小さい有界の誤差しか出さない:

摂動ケプラー(Kepler)問題

以下は摂動ケプラー問題として知られるハミルトン系をロードする.また,位置変数をハミルトン形式で定義して,積分区間と刻み幅を設定する:
この系には,H および L として定義される2つの不変量がある:

実験により,投影過程では使用できる不変量をすべて使うことが大切であることが分かる([HLW02]を参照).以下を使って得られた解を考えてみる.

ロトカ・ヴォルテラ系

ロトカ・ヴォルテラ系の制約射影の例は,「ロトカ・ヴォルテラ(LotkaVolterra)方程式の数値解法」に記載されている.

オイラー方程式

オイラー方程式の制約射影の例は,「剛体ソルバ」に記載されている.

オプションの要約

オプション名
デフォルト値
"Invariants"None微分方程式系の不変量を指定する
"IterationSafetyFactor"不変量の反復解法で使う安全率を指定する
MaxIterationsAutomatic不変量の反復解法の最大反復回数を指定する
Method"StiffnessSwitching"微分方程式系の数値積分に使うメソッドを指定する

"Projection"メソッドのオプション