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

MATLAB® で使用可能な質量行列ソルバーの 1 つを使用して微分代数方程式を解きます。このワークフローを使用するには、はじめに微分代数方程式 (DAE) の求解の手順 1、2、および 3 を完了してください。その後、ode15i の代わりに質量行列ソルバーを使用します。

この例では、ode15s または ode23t の使用方法を示します。他のソルバーの詳細については、ODE ソルバーの選択を参照してこのページのワークフローを適用してください。

手順 1. DAE の関数ハンドルへの変換

微分代数方程式 (DAE) の求解の手順 1 ～ 3 のコードは次のとおりです。

syms x(t) y(t) T(t) m r g
eqn1 = m*diff(x(t), 2) == T(t)/r*x(t);
eqn2 = m*diff(y(t), 2) == T(t)/r*y(t) - m*g;
eqn3 = x(t)^2 + y(t)^2 == r^2;
eqns = [eqn1 eqn2 eqn3];
vars = [x(t); y(t); T(t)];
origVars = length(vars);
[eqns,vars] = reduceDifferentialOrder(eqns,vars);
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);

reduceDAEIndex と reduceRedundancies の出力により、方程式のベクトル DAEs と変数のベクトル DAEvars が得られました。ode15s または ode23t を使用するには、2 つの関数ハンドルが必要です。1 つは DAE 系の質量行列を、もう 1 つは質量行列方程式の右辺を表します。これらの関数ハンドルは ODE 系の質量行列の等価な表現です。ここで $\mathit{M}\left(\mathit{t},\mathit{y}\left(\mathit{t}\right)\right){\mathit{y}}^{\prime }\left(\mathit{t}\right)=\mathit{f}\left(\mathit{t},\mathit{y}\left(\mathit{t}\right)\right)$ です。

massMatrixForm を使用してこれらの関数ハンドルを求め、質量行列 M と右辺 F を得ます。

[M,f] = massMatrixForm(DAEs,DAEvars)

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

f = 
$$\left(\begin{array}{c}\frac{T\left(t\right)\hspace{0.17em}x\left(t\right)-m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxtt}\left(t\right)}{r}\\ -\frac{g\hspace{0.17em}m\hspace{0.17em}r-T\left(t\right)\hspace{0.17em}y\left(t\right)+m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dytt}\left(t\right)}{r}\\ r^2-{x\left(t\right)}^2-{y\left(t\right)}^2\\ -2\hspace{0.17em}\mathrm{Dxt}\left(t\right)\hspace{0.17em}x\left(t\right)-2\hspace{0.17em}\mathrm{Dyt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ -2\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2-2\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2-2\hspace{0.17em}\mathrm{Dxtt}\left(t\right)\hspace{0.17em}x\left(t\right)-2\hspace{0.17em}\mathrm{Dytt}\left(t\right)\hspace{0.17em}y\left(t\right)\\ -\mathrm{Dytt}\left(t\right)\\ -\mathrm{Dyt}\left(t\right)\end{array}\right)$$

DAEs の方程式には変数のベクトル DAEvars で指定されていないシンボリック パラメーターが含まれても構いません。DAEs および DAEvars から得られた symvar の出力に対して setdiff を使用して、これらのパラメーターを求めます。

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars)

extraParams = $$\left(\begin{array}{ccc}g& m& r\end{array}\right)$$

質量行列 M にはこれらの余分なパラメーターはありません。したがって、odeFunction を使用して M を関数ハンドルに直接変換します。

M = odeFunction(M,DAEvars);

f を関数ハンドルに変換します。odeFunction への追加の入力としてパラメーターを別途指定します。

f = odeFunction(f,DAEvars,g,m,r);

以降のワークフローは純粋に数値的です。パラメーター値を設定して関数ハンドルを作成します。

gVal = 9.81;
mVal = 1;
rVal = 1;
F = @(t,Y) f(t,Y,gVal,mVal,rVal);

手順 2. 初期条件を求める

ソルバーは関数ハンドル内のあらゆる変数に対して初期値を必要とします。MATLAB 関数decicを使用して、方程式を充足する初期値を求めます。decic は推測された (方程式を満たさない可能性のある) 初期条件も取り扱うことができ、それらの推測を使用して充足する初期条件を求めようとします。decic は失敗することもあり、そのような場合は問題に対し整合性のある初期値を手動で設定する必要があります。

はじめに、DAEvars の変数をチェックします。

DAEvars

DAEvars = 
$$\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ T\left(t\right)\\ \mathrm{Dxt}\left(t\right)\\ \mathrm{Dyt}\left(t\right)\\ \mathrm{Dytt}\left(t\right)\\ \mathrm{Dxtt}\left(t\right)\end{array}\right)$$

ここで、Dxt(t) は x(t) の 1 次導関数、Dytt(t) は y(t) の 2 次導関数、以下同様です。7 × 1 ベクトルには 7 個の変数が含まれます。したがって、変数の初期値とその導関数の推測も 7 × 1 ベクトルでなければなりません。

振子の初期角変位が 30° つまり pi/6、座標の原点が振子のつり下げ点にあるとします。半径 rVal は 1 であるため、初期水平位置 x(t) は rVal*sin(pi/6) です。初期垂直位置 y(t) は -rVal*cos(pi/6) です。ベクトル y0est におけるこれらの変数の初期値を指定します。

残りの変数の初期値を任意に設定し、その導関数を 0 に設定します。これらは良い推測ではありません。しかしながら、この問題を充足します。問題において decic エラーが発生した場合は、より良い推測を与えてからdecicのページを参照してください。

y0est = [rVal*sin(pi/6); -rVal*cos(pi/6); 0; 0; 0; 0; 0];
yp0est = zeros(7,1);

質量行列 M と初期の推測 yp0est を含むオプションセットを作成し、数値検索の数値許容誤差を指定するオプション セットを指定します。

opt = odeset(Mass=M,InitialSlope=yp0est, ...
      RelTol=1e-7,AbsTol=1e-7);

MATLAB 関数 decic を使用して、変数とその導関数で整合性のある初期値を求めます。decic の最初の引数は、DAE を f(t,y,yp) = f(t,y,y') = 0 と記述した関数ハンドルでなければなりません。M および F の観点からすると、これは f(t,y,yp) = M(t,y)*yp - F(t,y) を意味します。

implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
[y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt)

y0 = 7×1

    0.4771
   -0.8788
   -8.6214
         0
    0.0000
   -2.2333
   -4.1135

yp0 = 7×1

         0
    0.0000
         0
         0
   -2.2333
         0
         0

ここで、系の質量行列 M および、導関数に対して整合性のある初期値のベクトル yp0 を含むオプション セットを作成します。系を解くときにこのオプションを使用します。

opt = odeset(opt,InitialSlope=yp0);

手順 3. DAE 系の求解

系を時間範囲 0 ≤ t ≤ 0.5 で積分して求解します。プロットにグリッド ラインと凡例を追加します。コードでは ode15s を使用します。代わりに、ode15s を ode23t に置き換えて ode23t を使用することもできます。

[tSol,ySol] = ode15s(F,[0, 0.5],y0,opt);
plot(tSol,ySol(:,1:origVars),"-o")

for k = 1:origVars
    S{k} = char(DAEvars(k));
end

legend(S,Location="best")
grid on

新しい変数を設定し、関数ハンドルと初期条件を再生成することで、異なるパラメーター値で系を解きます。

rVal に 2 を設定して、関数ハンドルと初期条件を再生成します。

rVal = 2;

F = @(t,Y) f(t,Y,gVal,mVal,rVal);
y0est = [rVal*sin(pi/6); -rVal*cos(pi/6); 0; 0; 0; 0; 0];
implicitDAE = @(t,y,yp) M(t,y)*yp - F(t,y);
[y0,yp0] = decic(implicitDAE,0,y0est,[],yp0est,[],opt);

opt = odeset(opt,InitialSlope=yp0);

新しいパラメーター値で系を解きます。

[tSol,ySol] = ode15s(F,[0, 0.5],y0,opt);
plot(tSol,ySol(:,1:origVars),"-o")

for k = 1:origVars
    S{k} = char(DAEvars(k));
end

legend(S,Location="best")
grid on

参考

トピック

• 微分代数方程式 (DAE) の求解
• 半線形 DAE 系の求解