半線形 DAE 系の求解

このワークフローは半線形微分代数方程式 (DAE) を解くための代替ワークフローで、標準ワークフローにおいて reduceDAEIndex が失敗して次の警告メッセージが発生したときだけ使用されます。The index of the reduced DAEs is larger than 1.標準ワークフローについては、微分代数方程式 (DAE) の求解を参照してください。

他の手順を始める前に微分代数方程式 (DAE) の求解の手順 1 と 2 を完了します。その後、手順 3 で reduceDAEIndex が失敗した場合、reduceDAEToODE を使用して微分指数を簡約します。reduceDAEToODE を使用する利点は、半線形 DAE を ODE (指数が 0 の DAE) に確実に簡約できることです。ただし、この関数は時間がかかり、半線形の DAE 系にのみ使用できます。系が半線形でない場合、reduceDAEToODE は失敗します。

DAE 系を解くにはこれらの手順を完了してください。

手順 1. 方程式と変数の指定

DAE 系は次のとおりです。

$$\begin{array}{l}m\frac{d^{2}x}{d{t}^2}=T\left(t\right)\frac{x\left(t\right)}{r}\\ m\frac{d^{2}y}{d{t}^2}=T\left(t\right)\frac{y\left(t\right)}{r}-mg\\ x{\left(t\right)}^2+y{\left(t\right)}^2=r^2\end{array}$$

syms を用いて独立変数と状態変数を指定します。

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);

手順 2. 微分の階数の簡約

DAE 系の "微分階数" とは、その方程式の最高微分階数です。MATLAB® を使用して DAE を解くには、微分階数を 1 に簡約する必要があります。ここで、1 番目および 2 番目の方程式には x(t) および y(t) の 2 次導関数が含まれています。したがって、微分階数は 2 です。

reduceDifferentialOrder を使用して系を 1 階の系に簡約します。関数 reduceDifferentialOrder は導関数を Dxt(t) や Dyt(t) のような新しい変数に置き換えます。eqns の式の右辺は 0 です。

[eqns,vars] = reduceDifferentialOrder(eqns,vars)

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

vars = 
$$\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)\end{array}\right)$$

手順 3. reduceDAEToODE を使った微分指数の簡約

eqns および vars により記述された DAE の微分指数を簡約するには、reduceDAEToODE を使用します。指数を簡約するため、reduceDAEToODE は新しい変数と方程式を系に追加します。reduceDAEToODE は制約も返します。これは、出力された ODE が初期の DAE と確実に等しい初期値を求めるのに役立つ条件です。

[ODEs,constraints] = reduceDAEToODE(eqns,vars)

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

constraints = 
$$\left(\begin{array}{c}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)\\ r^2-{x\left(t\right)}^2-{y\left(t\right)}^2\\ -2\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}{\mathrm{Dxt}\left(t\right)}^2-2\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}{\mathrm{Dyt}\left(t\right)}^2-2\hspace{0.17em}T\left(t\right)\hspace{0.17em}{x\left(t\right)}^2-2\hspace{0.17em}T\left(t\right)\hspace{0.17em}{y\left(t\right)}^2+2\hspace{0.17em}g\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}y\left(t\right)\end{array}\right)$$

手順 4. ODE から ode15s と ode23t の関数ハンドル

reduceDAEToODE の出力により、ODEs の方程式のベクトルと vars の変数のベクトルが得られます。ode15s または ode23t を使用するには、2 つの関数ハンドルが必要です。1 つは ODE 系の質量行列を、もう 1 つは質量行列方程式の右辺を含むベクトルを表します。これらの関数ハンドルは ODE 系の質量行列の等価な表現です。ここで、M(t,y(t))y’(t) = f(t,y(t)) です。

massMatrixForm を使用してこれらの関数ハンドルを求め、質量行列 massM (方程式の M) と右辺 f を得ます。

[massM,f] = massMatrixForm(ODEs,vars)

massM = 
$$\left(\begin{array}{ccccc}-1& 0& 0& 0& 0\\ 0& -1& 0& 0& 0\\ 0& 0& 0& m& 0\\ 0& 0& 0& 0& m\\ -4\hspace{0.17em}T\left(t\right)\hspace{0.17em}x\left(t\right)& 2\hspace{0.17em}g\hspace{0.17em}m\hspace{0.17em}r-4\hspace{0.17em}T\left(t\right)\hspace{0.17em}y\left(t\right)& -2\hspace{0.17em}{x\left(t\right)}^2-2\hspace{0.17em}{y\left(t\right)}^2& -4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dxt}\left(t\right)& -4\hspace{0.17em}m\hspace{0.17em}r\hspace{0.17em}\mathrm{Dyt}\left(t\right)\end{array}\right)$$

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

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

pODEs = symvar(ODEs);
pvars = symvar(vars);
extraParams = setdiff(pODEs,pvars)

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

別途指定する必要のあるパラメーターは質量 m、半径 r、および重力定数 g です。

odeFunction を使用して、massM および f を関数ハンドルに変換します。odeFunction への追加の入力として、シンボリック パラメーターを別途指定します。

massM = odeFunction(massM,vars,m,r,g);
f = odeFunction(f,vars,m,r,g);

以降のワークフローは純粋に数値的です。パラメーター値を設定して DAEs および constraints のパラメーター値を置き換えます。

mVal = 1;
rVal = 1;
gVal = 9.81;
ODEsNumeric = subs(ODEs,[m r g],[mVal rVal gVal]);
constraintsNumeric = subs(constraints,[m r g],[mVal rVal gVal]);

ode15s または ode23s の入力に適した関数ハンドルを作成します。

M = @(t,Y) massM(t,Y,mVal,rVal,gVal);
F = @(t,Y) f(t,Y,mVal,rVal,gVal);

手順 5. ode15s と ode23t の初期条件

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

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

vars

vars = 
$$\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)\end{array}\right)$$

ここで、Dxt(t) は x(t) の 1 次導関数です。5 × 1 ベクトルには 5 個の変数が含まれます。したがって、変数の初期値とその導関数の推測も 5 × 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];
yp0est = zeros(5,1);

系の質量行列 M を含み、数値検索の数値許容誤差を指定するオプション セットを作成します。

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

関数 decic を使用して ODE 系および代数制約に整合する初期値を求めます。この関数呼び出しのパラメーター [1,0,0,0,1] は y0est の最初と最後の要素を固定し、数値検索中に decic によって変更されないようにします。ここで、この修正は decic が充足するような初期条件を確実に求めるのに不可欠です。

[y0,yp0] = decic(ODEsNumeric,vars,constraintsNumeric,0, ...
                y0est,[1,0,0,0,1],yp0est,opt)

y0 = 5×1

    0.5000
   -0.8660
   -8.4957
   -0.0000
         0

yp0 = 5×1

   -0.0000
    0.0000
   -0.0000
   -4.2479
   -2.4525

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

opt = odeset(opt,InitialSlope=yp0);

手順 6. ode15s または ode23t による ODE 系の求解

系を時間範囲 0 ≤ t ≤ 0.5 で積分して求解します。プロットにグリッド ラインと凡例を追加します。ode15s を ode23s に置き換えてode23s を使用します。

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

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

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

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

rVal を 2 に設定し、手順を繰り返します。

rVal = 2;

ODEsNumeric = subs(ODEs,[m r g],[mVal rVal gVal]);
constraintsNumeric = subs(constraints,[m r g],[mVal rVal gVal]);
M = @(t,Y) massM(t,Y,mVal,rVal,gVal);
F = @(t,Y) f(t,Y,mVal,rVal,gVal);

y0est = [rVal*sin(pi/6); -rVal*cos(pi/6); 0; 0; 0];
opt = odeset(Mass=M,RelTol=1e-7,AbsTol=1e-7);
[y0, yp0] = decic(ODEsNumeric,vars,constraintsNumeric,0, ...
		 y0est,[1,0,0,0,1],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(vars(k));
end

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

参考

daeFunction | decic | findDecoupledBlocks | incidenceMatrix | isLowIndexDAE | massMatrixForm | odeFunction | reduceDAEIndex | reduceDAEToODE | reduceDifferentialOrder | reduceRedundancies

トピック

• 微分代数方程式 (DAE) の求解
• 質量行列ソルバーを使用した DAE の求解