Main Content

reduceDAEToODE

1 階半線形微分代数方程式系の同等の微分指数 0 の系への変換

ページ内をすべて折りたたむ

構文

newEqs = reduceDAEToODE(eqs,vars)
[newEqs,constraintEqs] = reduceDAEToODE(eqs,vars)
[newEqs,constraintEqs,oldIndex] = reduceDAEToODE(eqs,vars)

説明

newEqs = reduceDAEToODE(eqs,vars) は、高指数の 1 階半線形代数方程式系 eqs を同等の常微分方程式の系 newEqs に変換します。新しい系の微分指数は 0 です。つまり、vars の変数の導関数に対し newEqs のヤコビアンは可逆です。

例

[newEqs,constraintEqs] = reduceDAEToODE(eqs,vars) は拘束方程式のベクトルを返します。

例

[newEqs,constraintEqs,oldIndex] = reduceDAEToODE(eqs,vars) は元の半線形 DAE の系 eqs の微分指数 oldIndex を返します。

例

例

DAE 系の陰的 ODE 系への変換

微分代数方程式 (DAE) 系を陰的常微分方程式 (ODE) 系に変換します。

次の 2 つの微分代数方程式から成る系を作成します。ここで、シンボリック関数 x(t)、y(t) および z(t) は系の状態変数を表します。方程式と変数を 2 つのシンボリック ベクトル (方程式をシンボリック方程式のベクトル、変数をシンボリック関数呼び出しのベクトル) として指定します。

syms x(t) y(t) z(t)
eqs = [diff(x,t)+x*diff(y,t) == y,...
       x*diff(x, t)+x^2*diff(y) == sin(x),...
       x^2 + y^2 == t*z];
vars = [x(t), y(t), z(t)];

reduceDAEToODE を使用して微分指数が 0 になるように系を書き換えます。

newEqs = reduceDAEToODE(eqs, vars)
newEqs =
                            x(t)*diff(y(t), t) - y(t) + diff(x(t), t)
                diff(x(t), t)*(cos(x(t)) - y(t)) - x(t)*diff(y(t), t)
 z(t) - 2*x(t)*diff(x(t), t) - 2*y(t)*diff(y(t), t) + t*diff(z(t), t)

系を簡約して、簡約結果の詳細を返す

次の DAE 系の微分指数が低い (0 または 1) か高い (>1) か調べます。指数が 1 よりも高い場合、これを下げるためにまず reduceDAEIndex を、次に reduceDAEToODE を使用します。

微分代数方程式系を作成します。ここで、関数 x1(t)、x2(t) および x3(t) は系の状態変数を表します。また、系には関数 q1(t)、q2(t) および q3(t) が含まれます。これらの関数は状態関数を表しません。方程式と変数を 2 つのシンボリック ベクトル (方程式をシンボリック方程式のベクトル、変数をシンボリック関数呼び出しのベクトル) として指定します。

syms x1(t) x2(t) x3(t) q1(t) q2(t) q3(t)
eqs = [diff(x2) == q1 - x1,
       diff(x3) == q2 - 2*x2 - t*(q1-x1),
       q3 - t*x2 - x3];
vars = [x1(t), x2(t), x3(t)];

isLowIndexDAE を使用して系の微分指数をチェックします。この系に対して isLowIndexDAE は 0 (false) を返します。これは系の微分指数が 2 以上であることを意味します。

isLowIndexDAE(eqs, vars)
ans =
  logical
     0

系をはじめて書き換えようとする場合は、reduceDAEIndex を使用し、微分指数が 1 になるようにします。この系では、reduceDAEIndex が系の微分指数を 0 または 1 まで減らせないため、警告が出されます。

[newEqs, newVars] = reduceDAEIndex(eqs, vars)
Warning: Index of reduced DAEs is larger than 1.
 
newEqs =
                      x1(t) - q1(t) + diff(x2(t), t)
       Dx3t(t) - q2(t) + 2*x2(t) + t*(q1(t) - x1(t))
                             q3(t) - x3(t) - t*x2(t)
 diff(q3(t), t) - x2(t) - t*diff(x2(t), t) - Dx3t(t)
 
newVars =
   x1(t)
   x2(t)
   x3(t)
 Dx3t(t)

reduceDAEIndex で半線形系の指数を 0 または 1 に簡約できない場合は、reduceDAEToODE を試します。この関数の処理は大幅に遅い場合があるため、最初の選択としてはお勧めできません。2 つの出力引数をもつ構文の使用により拘束方程式も返されます。

[newEqs, constraintEqs] = reduceDAEToODE(eqs, vars)
newEqs =
                                             x1(t) - q1(t) + diff(x2(t), t)
                       2*x2(t) - q2(t) + t*q1(t) - t*x1(t) + diff(x3(t), t)
 diff(x1(t), t) - diff(q1(t), t) + diff(q2(t), t, t) - diff(q3(t), t, t, t)

constraintEqs =
 x1(t) - q1(t) + diff(q2(t), t) - diff(q3(t), t, t)
                            x3(t) - q3(t) + t*x2(t)
                     x2(t) - q2(t) + diff(q3(t), t)

3 つの出力引数をもつ構文を使用して、新しい方程式、拘束方程式および元の系 eqs の微分指数を返します。

[newEqs, constraintEqs, oldIndex] = reduceDAEToODE(eqs, vars)
newEqs =
                                             x1(t) - q1(t) + diff(x2(t), t)
                       2*x2(t) - q2(t) + t*q1(t) - t*x1(t) + diff(x3(t), t)
 diff(x1(t), t) - diff(q1(t), t) + diff(q2(t), t, t) - diff(q3(t), t, t, t)

constraintEqs =
 x1(t) - q1(t) + diff(q2(t), t) - diff(q3(t), t, t)
                            x3(t) - q3(t) + t*x2(t)
                     x2(t) - q2(t) + diff(q3(t), t)

oldIndex =
     3

入力引数

すべて折りたたむ

半線形 1 階 DAE の系。シンボリック方程式またはシンボリック式のベクトルとして指定します。

状態変数。x(t) など、シンボリックな関数または関数呼び出しのベクトルとして指定します。

例: [x(t),y(t)] または [x(t);y(t)]

出力引数

すべて折りたたむ

陰的常微分方程式系。シンボリック式の列ベクトルとして返されます。この系の微分指数は 0 です。

系の簡約化時に発生する拘束方程式。シンボリック式の列ベクトルとして返されます。これらの式は変数 vars に依存しますが、その導関数には依存しません。制約は newEqs の微分方程式の保存量です。つまり各制約の時間微分は newEqs の方程式の剰余を消滅させます。

これらの方程式を使用して DAE 系に整合する初期条件を決定できます。

元の DAE 系 eqs の微分指数。整数として返されます。

アルゴリズム

reduceDAEToODE の実装はガウスの消去法に基づいています。このアルゴリズムは reduceDAEIndex で使用される Pantelides アルゴリズムよりもより信頼性がありますが、大幅に時間がかかります。

バージョン履歴

R2014b で導入

参考

| | | | | | | | |

トピック