massMatrixForm

質量行列と、半線形連立微分代数方程式の右辺の抽出

構文

[M,F] = massMatrixForm(eqs,vars)

説明

[M,F] = massMatrixForm(eqs,vars) は、質量行列 M および 1 階半線形微分代数方程式 (DAE) 系の方程式の右辺 F を返します。vars の変数の導関数を一切含まない eqs の代数方程式はその質量行列 M の空の行に対応します。

質量行列 M および方程式の右辺 F はこの形で表されます。

$$M\left(t,x\left(t\right)\right)\dot{x}\left(t\right)=F\left(t,x\left(t\right)\right).$$

例

DAE 系の質量行列形式への変換

半線形連立微分代数方程式を質量行列形式に変換します。

次の連立微分代数方程式を作成します。ここで、関数 x(t) および y(t) は方程式の状態変数を表します。また、方程式にはシンボリック パラメーター alpha と beta、および関数 f(t,x,y) も含まれます。方程式と変数を 2 つのシンボリック ベクトル (方程式をシンボリック方程式のベクトル、変数をシンボリック関数呼び出しのベクトル) として指定します。

syms x(t) y(t) f(t,x,y) alpha beta
eqs = [x(t) - diff(x(t),t) == alpha*f(t,x(t),y(t)); ...
       y(t) + diff(y(t),t) == beta*f(t,x(t),y(t))]

eqs = 
$$\left(\begin{array}{c}x\left(t\right)-\frac{\partial }{\partial t}\mathrm{ }x\left(t\right)=\alpha \hspace{0.17em}f\left(t,x\left(t\right),y\left(t\right)\right)\\ \frac{\partial }{\partial t}\mathrm{ }y\left(t\right)+y\left(t\right)=\beta \hspace{0.17em}f\left(t,x\left(t\right),y\left(t\right)\right)\end{array}\right)$$

vars = [x(t),y(t)];

この方程式の質量行列形式を求めます。

[M,F] = massMatrixForm(eqs,vars)

M = 
$$\left(\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right)$$

F = 
$$\left(\begin{array}{c}\alpha \hspace{0.17em}f\left(t,x\left(t\right),y\left(t\right)\right)-x\left(t\right)\\ \beta \hspace{0.17em}f\left(t,x\left(t\right),y\left(t\right)\right)-y\left(t\right)\end{array}\right)$$

数値ソルバー ode15s を使用してこの方程式を求解します。ode15s を使用する前に、方程式のシンボリック パラメーターに値 alpha = 0.01、beta = 0.02、f(t,x,y) = x*y を代入します。また、状態変数 x(t) および y(t) を matlabFunction に受け入れ可能な変数 Y1 および Y2 に置き換えます。

syms Y1 Y2
M = subs(M,[vars,alpha,beta,f],[Y1,Y2,0.01,0.02,@(t,x,y) x*y]);
F = subs(F,[vars,alpha,beta,f],[Y1,Y2,0.01,0.02,@(t,x,y) x*y]);

次の関数ハンドル MM および FF を作成します。これらの関数ハンドルを odeset および ode15s の入力引数として使用できます。これらの関数では状態変数が列ベクトルとして指定されている必要があります。

MM = matlabFunction(M,Vars={t,[Y1; Y2]});
FF = matlabFunction(F,Vars={t,[Y1; Y2]});

ode15s を使用して連立方程式を解きます。変数 Y1 および Y2 の初期値を 20 および 10 として指定します。時間間隔 [0 10] にわたり、解をプロットします。

opt = odeset(Mass=MM);
ode15s(FF,[0 15],[50; 20],opt)

状態変数のプログラムによる特定と生成

次の連立微分代数方程式を作成します。ここで、シンボリック関数 x1(t)、x2(t) および x3(t) は方程式の状態変数を表します。また、方程式にはシンボリック パラメーター a および b が含まれます。

syms x1(t) x2(t) x3(t) a b
eqs = [a*diff(x1,t) == x1*x3- x2;
       diff(x2,t) == b*x1 - 1;
       0 == x1 + x2 + x3]

eqs(t) = 
$$\left(\begin{array}{c}a\hspace{0.17em}\frac{\partial }{\partial t}\mathrm{ }{x}_1\left(t\right)=x_1\left(t\right)\hspace{0.17em}{x}_3\left(t\right)-x_2\left(t\right)\\ \frac{\partial }{\partial t}\mathrm{ }{x}_2\left(t\right)=b\hspace{0.17em}{x}_1\left(t\right)-1\\ 0=x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\end{array}\right)$$

連立方程式の状態変数を求めます。多くの場合、これらの変数は時間のみの関数です。状態変数をプログラムによって特定するには、symfunOf を t として指定して findSymType を使用します。

vars = findSymType(eqs,"symfunOf",t)

vars = $$\left(\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)\end{array}\right)$$

この方程式の質量行列形式を求めます。

[M,F] = massMatrixForm(eqs,vars)

M = 
$$\left(\begin{array}{ccc}a& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)$$

F = 
$$\left(\begin{array}{c}{x}_1\left(t\right)\hspace{0.17em}{x}_3\left(t\right)-x_2\left(t\right)\\ b\hspace{0.17em}{x}_1\left(t\right)-1\\ x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\end{array}\right)$$

次に、a および b の値をそれぞれ 2 および 1 として代入します。また、状態変数を matlabFunction に受け入れ可能な変数 Y1、Y2 などに置き換えます。sym を使用して変数 Y1、Y2 などをプログラムによって生成できます。

numVars = numel(vars);
vecY = sym("Y",[1,numVars]);
M = subs(M,[a b vars],[2 1 vecY]);
F = subs(F,[a b vars],[2 1 vecY]);

次の関数ハンドル MM および FF を作成します。これらの関数ハンドルを odeset および ode15s の入力引数として使用します。これらの関数では状態変数が列ベクトルとして指定されている必要があります。vecY のシンボリック変数は複素数であるため、非共役転置を使用して vecY を列ベクトルに変換します。

MM = matlabFunction(M,Vars={t,vecY.'});
FF = matlabFunction(F,Vars={t,vecY.'});

ode15s を使用して連立方程式を解きます。変数 Y1、Y2、および Y3 の初期値をそれぞれ 1、1、および -2 として指定します。時間間隔 [0 15] にわたり、解をプロットします。

opt = odeset(Mass=MM);
ode15s(FF,[0 15],[1; 1; -2],opt)

入力引数

eqs — 半線形 1 階 DAE 系

シンボリック方程式のベクトル | シンボリック式のベクトル

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

vars — 状態変数

シンボリック関数のベクトル | シンボリック関数呼び出しのベクトル

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

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

出力引数

M — 質量行列

シンボリック行列

方程式の質量行列。シンボリック行列として返されます。行数は eqs の方程式の数に対応し、列数は vars の変数の数に対応します。

F — 方程式の右辺

シンボリック式のシンボリック列ベクトル

方程式の右辺。シンボリック式の列ベクトルとして返されます。このベクトルの要素の数は eqs の方程式の数と一致します。

バージョン履歴

R2014b で導入

参考

findDecoupledBlocks | daeFunction | decic | incidenceMatrix | isLowIndexDAE | matlabFunction | odeFunction | ode15s | odeset | reduceDAEIndex | reduceDAEToODE | reduceDifferentialOrder | reduceRedundancies

トピック

• 質量行列ソルバーを使用した DAE の求解