Main Content

findDecoupledBlocks

方程式系の分離ブロックの検索

説明

[eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs,vars) は、変数のサブセットの定義に使用可能な方程式のサブセット (ブロック) を特定します。変数の数 vars は方程式の数 eqs と一致しなければなりません。

i 番目のブロックは、vars(varsBlocks{i}) の変数を決定する一連の方程式です。vars([varsBlocks{1},…,varsBlocks{i-1}]) の変数は前の方程式のブロックによって再帰的に決定されます。1 番目の方程式のブロックを解いて変数の 1 番目のブロックを求めた後、eqs(eqsBlocks{2}) によって与えられた方程式の 2 番目のブロックは、変数の 2 番目のブロック (vars(varsBlock{2})) によって与えられた変数のサブセットおよび 1 番目のブロックの変数 (これらの変数はこの時点で判明している) のみを含む、分離された方程式のサブセットを定義します。したがって、非自明なブロック分解が可能な場合、多くの変数を含む大きな方程式系の求解プロセスを、いくつかのより小さなサブの系のステップに分割することができます。

ブロック数 length(eqsBlocks)length(varsBlocks) と一致します。length(eqsBlocks) = length(varsBlocks) = 1 の場合、方程式の非自明なブロック分解は不可能です。

すべて折りたたむ

シンボリックの微分代数方程式 (DAE) 系のブロック下三角行列 (BLT) 分解を計算します。

次の 4 つの微分代数方程式から成る系を作成します。ここで、シンボリック関数呼び出し x1(t)x2(t)x3(t) および x4(t) は、系の状態変数を表します。また、系にはシンボリック パラメーター c1c2c3c4、および関数 f(t,x,y)g(t,x,y) が含まれます。

syms x1(t) x2(t) x3(t) x4(t)
syms c1 c2 c3 c4
syms f(t,x,y) g(t,x,y)

eqs = [c1*diff(x1(t),t)+c2*diff(x3(t),t)==c3*f(t,x1(t),x3(t));...
       c2*diff(x1(t),t)+c1*diff(x3(t),t)==c4*g(t,x3(t),x4(t));...
       x1(t)==g(t,x1(t),x3(t));...
       x2(t)==f(t,x3(t),x4(t))];

vars = [x1(t),x2(t),x3(t),x4(t)];

findDecoupledBlocks を使用して、系のブロック構造体を見つけます。

[eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs,vars)
eqsBlocks=1×3 cell array
    {[1 3]}    {[2]}    {[4]}

varsBlocks=1×3 cell array
    {[1 3]}    {[4]}    {[2]}

最初のブロックには 2 つの変数の 2 つの方程式が含まれます。

eqs(eqsBlocks{1})
ans = 

(c1t x1(t)+c2t x3(t)=c3f(t,x1(t),x3(t))x1(t)=g(t,x1(t),x3(t)))

vars(varsBlocks{1})
ans = (x1(t)x3(t))

このブロックの変数 x1(t)x3(t) の解を得た後、次の方程式ブロックを求解できます。このブロックは 1 つの方程式で構成されています。

eqs(eqsBlocks{2})
ans = 

c2t x1(t)+c1t x3(t)=c4g(t,x3(t),x4(t))

このブロックには 1 つの変数が関与しています。

vars(varsBlocks{2})
ans = x4(t)

ブロック 2 の方程式の変数 x4(t) の解を得た後、残りの方程式ブロック eqs(eqsBlocks{3}) で残りの変数 vars(varsBlocks{3}) を定義します。

eqs(eqsBlocks{3})
ans = x2(t)=f(t,x3(t),x4(t))
vars(varsBlocks{3})
ans = x2(t)

系をブロック下三角行列形式に変換する順列を求めます。

eqsPerm = [eqsBlocks{:}]
eqsPerm = 1×4

     1     3     2     4

varsPerm = [varsBlocks{:}]
varsPerm = 1×4

     1     3     4     2

系を方程式のブロック下三角行列に変換します。

eqs = eqs(eqsPerm)
eqs = 

(c1t x1(t)+c2t x3(t)=c3f(t,x1(t),x3(t))x1(t)=g(t,x1(t),x3(t))c2t x1(t)+c1t x3(t)=c4g(t,x3(t),x4(t))x2(t)=f(t,x3(t),x4(t)))

vars = vars(varsPerm)
vars = (x1(t)x3(t)x4(t)x2(t))

生成された系の接続行列を求めます。接続行列は、並べ替えられた方程式系にサイズが 22 列、11 列、11 列の 3 つの対角ブロックがあることを示しています。

incidenceMatrix(eqs, vars)
ans = 4×4

     1     1     0     0
     1     1     0     0
     1     1     1     0
     0     1     1     1

線形代数系で方程式のブロックを求め、方程式の各ブロックを 1 番目のものから順に解くことにより系を求解します。

次の線形代数方程式系を作成します。

syms x1 x2 x3 x4 x5 x6 c1 c2 c3

eqs = [c1*x1 + x3 + x5 == c1 + c2 + 1;...
       x1 + x3 + x4 + 2*x6 == 4 + c2;...
       x1 + 2*x3 + c3*x5 == 1 + 2*c2 + c3;...
       x2 + x3 + x4 + x5 == 2 + c2;...
       x1 - c2*x3 + x5 == 2 - c2^2;...
       x1 - x3 + x4 - x6 == 1 - c2];

vars = [x1,x2,x3,x4,x5,x6];

findDecoupledBlocks を使用して系を下三角行列形式に変換します。この系では、findDecoupledBlocks は 3 つの方程式ブロックと対応する変数を識別します。

[eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs,vars)
eqsBlocks=1×3 cell array
    {[1 3 5]}    {[2 6]}    {[4]}

varsBlocks=1×3 cell array
    {[1 3 5]}    {[4 6]}    {[2]}

1 番目のブロックの変数を識別します。このブロックは 3 つの変数の 3 つの方程式で構成されています。

vars(varsBlocks{1})
ans = (x1x3x5)

方程式の 1 番目のブロックを変数の 1 番目のブロックに対して求解し、解を対応する変数に代入します。

[x1,x3,x5] = solve(eqs(eqsBlocks{1}),vars(varsBlocks{1}))
x1 = 1
x3 = c2
x5 = 1

2 番目のブロックの変数を識別します。このブロックは 2 つの変数の 2 つの方程式から構成されています。

vars(varsBlocks{2})
ans = (x4x6)

この方程式のブロックを解き、解を対応する変数に代入します。

[x4,x6] = solve(eqs(eqsBlocks{2}),vars(varsBlocks{2}))
x4 = 

x33-x1-c23+2

x6 = 

2c23-2x33+1

subs を使用して変数 x1x3、および x5 の既知の値に対する結果を評価します。

x4 = subs(x4)
x4 = 1
x6 = subs(x6)
x6 = 1

3 番目のブロックの変数を識別します。このブロックは 1 つの変数の 1 つの方程式で構成されています。

vars(varsBlocks{3})
ans = x2

x2 に解を代入してこの方程式を解きます。

x2 = solve(eqs(eqsBlocks{3}),vars(varsBlocks{3}))
x2 = c2-x3-x4-x5+2

subs を使用して、この系の他のすべての変数の既知の値に対する結果を評価します。

x2 = subs(x2)
x2 = 0

あるいは、for ループを使用してこの例を書き換えることもできます。これにより、この例をより大きな方程式系に使用することができます。

syms x1 x2 x3 x4 x5 x6 c1 c2 c3

eqs = [c1*x1 + x3 + x5 == c1 + c2 + 1;...
       x1 + x3 + x4 + 2*x6 == 4 + c2;...
       x1 + 2*x3 + c3*x5 == 1 + 2*c2 + c3;...
       x2 + x3 + x4 + x5 == 2 + c2;...
       x1 - c2*x3 + x5 == 2 - c2^2
       x1 - x3 + x4 - x6 == 1 - c2];

vars = [x1,x2,x3,x4,x5,x6];

[eqsBlocks,varsBlocks] = findDecoupledBlocks(eqs,vars);

vars_sol = vars;

for i = 1:numel(eqsBlocks)
  sol = solve(eqs(eqsBlocks{i}),vars(varsBlocks{i}));
  vars_sol_per_block = subs(vars(varsBlocks{i}),sol);
  for k=1:i-1
    vars_sol_per_block = subs(vars_sol_per_block,vars(varsBlocks{k}),...
                         vars_sol(varsBlocks{k}));
  end
  vars_sol(varsBlocks{i}) = vars_sol_per_block
end
vars_sol = (1x2c2x41x6)
vars_sol = (1x2c2111)
vars_sol = (10c2111)

入力引数

すべて折りたたむ

方程式系。シンボリックな方程式または式のベクトルとして指定します。

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

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

出力引数

すべて折りたたむ

方程式のブロックを定義するインデックス。cell 配列として返されます。各ブロックのインデックスは倍精度整数の行ベクトルです。方程式の i 番目のブロックは方程式 eqs(eqsBlocks{i}) で構成され、vars(varsBlocks{1:i}) の変数のみを含みます。

変数のブロックを定義するインデックス。cell 配列として返されます。各ブロックのインデックスは倍精度整数の行ベクトルです。方程式の i 番目のブロックは方程式 eqs(eqsBlocks{i}) で構成され、vars(varsBlocks{1:i}) の変数のみを含みます。

ヒント

  • 実装されているアルゴリズムでは、vars の各変数について、eqs にはその変数を含む方程式が、少なくとも 1 つは必ずあることが要求されます。同じ方程式が他の変数にも一致するようなことがあってはなりません。系がこの条件を満たさない場合、findDecoupledBlocks はエラーをスローします。特に、findDecoupledBlockslength(eqs) = length(vars) であることを要求します。

  • 順列 e = [eqsBlocks{:}] をベクトル eqs に、v = [varsBlocks{:}] をベクトル vars に適用すると、ブロック下三角行列スパース パターンをもつ接続行列 incidenceMatrix(eqs(e), vars(v)) が生成されます。

バージョン履歴

R2014b で導入