微分代数方程式系の状態空間法

ブロック下三角(BLT)形式オプションのまとめ

次の形式の分割された微分代数方程式系を考えてみる．

ヤコビ行列 および が両方可逆ならば，陰関数定理により，その系は状態空間形式で以下のように表すことができる．

可逆条件は実質的に系が指数1であることと同等なので，状態空間形式に基づくメソッドは指数1の微分代数方程式系に適している．最も一般的な形式は であり， が恒等行列であるものである．

時間積分メソッドの"StateSpace"はこの表現に基づいている． の値が与えられると，および を見付けるためにニュートンの反復法が使われる．右辺の関数の評価だけを必要する，根底の常微分方程式メソッドに対して，値 が返される．評価と評価の間に，次の評価の反復を初期化するために，前に計算された導関数の値と代数的変数が保存される．

状態空間形式から得られる常微分方程式は，非常に硬いことが多いためヤコビアンを必要とする．これは有限差分を用いて計算することができるが，行列 と がすでに計算されているために，問題構造を使った方がより正確で効率的である．右辺のヤコビアンは についての分割された系を微分することにより計算することができる．

右辺のヤコビアンを解くと，下が得られる．

および の行列分解は，通常ニュートン反復から保存されるので，転置の適用は効率的に計算することができる．

"StateSpace"メソッドは，各時間刻みに対してではなく各評価に対して，時によっては繰返し と を解かなければならず，1つの時間刻みに対して多数の評価があるかもしれないので，通常デフォルトの"IDA"マルチステップメソッドよりも遅い． 場合によってはスピードをオフセットする可能性のある利点も存在する．その利点とは，関数評価に対して局所的な反復はより安定しているということ，A安定の内在する積分法は全体の安定性を向上させるために使うことができるということ，内在するメソッドが1ステップであればそのメソッドは実質的に1ステップ法であるということ，メソッドは任意精度の解を許可するために実装されるということである．

"StateSpace"メソッドのMethodオプションを使うと，内在する積分メソッドを選ぶことができる．

ブロック下三角(BLT)形式

系の構造によっては，その系をブロック下三角形式に並べ替えた方がよい場合がある．ブロック下三角形式を使うことで，小さいサイズで数の多い部分問題の解となる．最大のブロックサイズを大きく減少させることができるならば，この構造を使うことでニュートン反復法の計算の複雑性を大きく減少させることができる．

ブロック下三角形式は，矛盾のない初期化を含め，微分代数方程式系を解く方法，指数減少法のダミー導関数を設定する方法，疎な線形系を解く方法の多くで使われている．ここでの説明は，状態空間法での使用に焦点を当てるが，他にも当てはまることもある．

ブロック下三角形式の最も簡単な形は，行列アルゴリズムであり，それがWolfram言語でアクセスするときの最も一般的な方法である．

SparseArray`BlockTriangularOrdering[s] s[[rk,ck]]が並べ替えられた行列s[[r,c]]の対角に沿った k 番目のブロックとなるように，行列 s の行と列の並べ替え{{r1,…,rb},{c1,…,cb}}を与える．ここで r=Flatten[{r1,…,rb}]，c=Flatten[{c1,…,cb}]，並べ替えられた行列の他のすべての要素は対角要素より下に現れる

ブロック下三角順を得る．

s が特異でない場合，対角要素を上から見ていくと，それぞれの連続するブロックの変数は前のブロックの変数にのみ依存する．つまり，線形方程式 は対角を上から順に，各ブロック部分系を順に解くことで解くことができ，計算努力を大きく節約することができる場合がある．

これを，小さい例で見てみよう．

次は，パイプで連結されたタンクの高さhと圧力pの変化を記述するNDSolveの例題からの線形微分代数方程式系である：
{h′,p}についてヤコビアンを得る：
ブロック下三角形式にする：
並べ替えられた行列を表示する：
並べ替えられた変数と方程式を表示する：

この系の場合，ブロックはすべてサイズ1であるため，特に解きやすい．ブロック下三角順の方程式のリストから分かるように，変数値を連続して選べるからである．状態空間メソッドの場合，h={h1,h2,h3}が与えられるので導関数h′を得ることは自明である．

ブロック下三角行列の便利さは，非線形問題にも当然当てはまる．

以下は，制約条件付きの振り子問題のダミー導関数を持つ指数減少法に由来する，指数1の微分代数方程式系である：
方程式に現れる最高階の導関数についてヤコビアンを計算する：
ブロック下三角順を得る：
並べ替えられた行列を表示する：

この形式から，ならばヤコビアンは特異であることが明白である．

並べ替えられた変数と方程式を表示する：

前と同様に，対角線を上から下に向かって，部分問題を順に解いていくことができる．x'，y，yd1が見付かると（非線形問題は，スカラーのニュートンソルバで実行することができる），x1'，yd2， λ を得るためには3×3のブロックだけを解けばよい．

方程式 には に対する解が2つある．ニュートン反復法を開始するために，前の値を使うことにより連続性が維持される．

メソッドオプション"BlockLowerTriangular"を使うと，さまざまな方法でブロック下三角順を使うよう状態空間メソッドに指示することができる．

Automaticどのようにブロック下三角順を使うかを自動的に決定する
Falseブロック下三角順を使わない
"Ordering"ブロック下三角順を使って，変数と方程式を並べる
"Solve"連続した部分系を解くためにブロック下三角順を使うことで解く

"BlockLowerTriangular"メソッドオプションの設定

でヤコビアンの特異値を避ける初期条件：
ブロック下三角順を使わずに状態空間メソッドを使って多数の区間上で解く：
ブロック下三角順で方程式と変数を並べる：
評価では，ブロック下三角の対角成分に従った部分系により解く：
解をプロットする：

この例では，系が大変小さいため，時間が非常に近い．ブロック下三角順のヤコビアンの小さいブロックを持った大きい系では，その差分がもっとはっきりしている．

オプションのまとめ

オプション名
デフォルト値
"BlockLowerTriangular"Automaticブロック下三角順の使用を指定する
MethodAutomatic内在する常微分方程式の時間積分法を指定する

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