ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

微分代数方程式 (DAE) の求解

MATLAB® と Symbolic Math Toolbox™ を使用して微分代数方程式 (DAE) の解を求めることができます。

微分代数方程式 (関数、あるいは状態変数、x(t)=[x1(t),,xn(t)] を含む) は、次のような形をもちます

F(t,x(t),x˙(t))=0

ここで、t は独立変数です。方程式の数 F=[F1,,Fn] は状態変数の数 x(t)=[x1(t),,xn(t)] と一致しなければなりません。

ほとんどの DAE 系は ode15i のような MATLAB ソルバーの入力として直接与えるのに適した形式になっていないため、Symbolic Math Toolbox の機能を用いてはじめに適切な形式に変換します。この機能は DAE の微分インデックス (系を ODE に換算するのに必要な微分の数) を 1 または 0 に減らしてから、DAE 系を MATLAB ソルバーに適した数値関数ハンドルに変換します。今度は、MATLAB ソルバーである ode15iode15s または ode23t のいずれかを使用して DAE を求解します。

これらの手順を完了することで DAE 系を解きます。

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

ここでは振子の DAE を解くことで DAE ワークフローを示します。

状態変数は次のとおりです。

  • 振子の水平位置 x(t)

  • 振子の垂直位置 y(t)

  • 振子の飛び出しを妨げる力 T(t)

変数は次のとおりです。

  • 振子の質量 m

  • 振子の長さ r

  • 重力定数 g

DAE 系は次のとおりです。

md2xdt2=T(t)x(t)rmd2ydt2=T(t)y(t)rmgx(t)2+y(t)2=r2

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. 微分の階数の縮小

2.1 (オプション) 変数の発生範囲の確認

この手順は "オプション" です。接続行列を見ることで DAE 系での変数の発生箇所を確認できます。この手順により、入力からは発生しないすべての変数を発見し、vars ベクトルから取り除くことができます。

incidenceMatrix を使用して新しい系の接続行列を表示します。incidenceMatrix の出力は各方程式を行に、各変数を列にもちます。この系には 3 つの方程式と 3 つの状態変数が含まれるため、incidenceMatrix33 列の行列を返します。行列には 10 があります。ここで 1 は状態変数の出現を表します。たとえば、位置 (2,3)1 は、2 番目の方程式には 3 番目の状態変数 T(t) が含まれることを意味します。

M = incidenceMatrix(eqns, vars)
M =
     1     0     1
     0     1     1
     1     1     0

接続行列のある列がすべて 0 だった場合は、その状態変数が DAE 系に出現しないため取り除かなければなりません。

2.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 =
       m*diff(Dxt(t), t) - (T(t)*x(t))/r
 g*m + m*diff(Dyt(t), t) - (T(t)*y(t))/r
                   x(t)^2 + y(t)^2 - r^2
                  Dxt(t) - diff(x(t), t)
                  Dyt(t) - diff(y(t), t)

vars =
   x(t)
   y(t)
   T(t)
 Dxt(t)
 Dyt(t)

オプションで、導関数が新しい変数にどのように対応するかを示すには、reduceDifferentialOrder を 3 つの出力引数で呼び出します。系の新しい変数の場所を見つけるには、incidenceMatrix を使用します。

手順 3. 微分指数の確認と縮小

3.1 系の微分指数の確認

isLowIndexDAE を使用して DAE 系の微分指数を確認します。指数が 0 または 1 である場合、isLowIndexDAE は論理値 1 (true) を返し、手順 3.2 を飛ばして手順 4. DAE 系の MATLAB 関数ハンドルへの変換に進むことができます。ここで、isLowIndexDAE が論理値 0 (false) を返します。これは微分指数が 1 よりも大きいため、縮小する必要があることを意味します。

isLowIndexDAE(eqns,vars)
ans =
  logical
     0

3.2 reduceDAEIndex による微分指数の縮小

微分指数を縮小するために、関数 reduceDAEIndex は与えられた方程式から導出される新しい方程式を追加し、高次の導関数を新しい変数に置き換えます。reduceDAEIndex が失敗して警告が発生した場合は、ワークフロー半線形 DAE 系の求解に記載されているとおり、代替となる関数 reduceDAEToODE を使用します。

eqns および vars の DAE の微分指数を縮小します。

[DAEs,DAEvars] = reduceDAEIndex(eqns,vars)
DAEs =
                                             m*Dxtt(t) - (T(t)*x(t))/r
                                       g*m + m*Dytt(t) - (T(t)*y(t))/r
                                                 x(t)^2 + y(t)^2 - r^2
                                                      Dxt(t) - Dxt1(t)
                                                      Dyt(t) - Dyt1(t)
                                       2*Dxt1(t)*x(t) + 2*Dyt1(t)*y(t)
 2*Dxt1t(t)*x(t) + 2*Dxt1(t)^2 + 2*Dyt1(t)^2 + 2*y(t)*diff(Dyt1(t), t)
                                                    Dxtt(t) - Dxt1t(t)
                                            Dytt(t) - diff(Dyt1(t), t)
                                               Dyt1(t) - diff(y(t), t)

DAEvars =
     x(t)
     y(t)
     T(t)
   Dxt(t)
   Dyt(t)
  Dytt(t)
  Dxtt(t)
  Dxt1(t)
  Dyt1(t)
 Dxt1t(t)

メモ

reduceDAEIndex が失敗して警告が発生した場合は、半線形 DAE 系の求解に記述されたワークフローを代わりに使用します。

多くの場合、reduceDAEIndex は削除できる冗長な方程式や変数を挿入します。reduceRedundancies を使用して冗長な方程式と変数を削除します。

[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)
DAEs =
                              -(T(t)*x(t) - m*r*Dxtt(t))/r
                       (g*m*r - T(t)*y(t) + m*r*Dytt(t))/r
                                     x(t)^2 + y(t)^2 - r^2
                             2*Dxt(t)*x(t) + 2*Dyt(t)*y(t)
 2*Dxtt(t)*x(t) + 2*Dytt(t)*y(t) + 2*Dxt(t)^2 + 2*Dyt(t)^2
                                 Dytt(t) - diff(Dyt(t), t)
                                    Dyt(t) - diff(y(t), t)
DAEvars =
    x(t)
    y(t)
    T(t)
  Dxt(t)
  Dyt(t)
 Dytt(t)
 Dxtt(t)

新しい系の微分指数を確認します。これで isLowIndexDAE は論理値 1 (true) を返すようになります。つまり、系の微分指数は 0 または 1 です。

isLowIndexDAE(DAEs,DAEvars)
ans =
  logical
     1

手順 4. DAE 系の MATLAB 関数ハンドルへの変換

この手順では MATLAB ODE ソルバーである ode15i 用の関数ハンドルを作成します。これは汎用ソルバーです。ode15sode23t のような質量行列に特化したソルバーを使用するには、質量行列ソルバーを使用した DAE の求解を参照してください。ODE ソルバーの選択 (MATLAB)も参照してください。

reduceDAEIndexDAEs の方程式のベクトルと DAEvars の変数のベクトルを出力します。ode15i を使用するには、DAE 系を表す関数ハンドルが必要です。

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

pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars)
extraParams =
[ g, m, r]

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

daeFunction を使用して関数ハンドルを作成します。daeFunction への追加の入力引数として、シンボリック パラメーターを別途指定します。

f = daeFunction(DAEs, DAEvars, g, m, r);

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

g = 9.81;
m = 1;
r = 1;
F = @(t, Y, YP) f(t, Y, YP, g, m, r);

手順 5. ソルバーの初期条件を求める

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

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

DAEvars
DAEvars =
    x(t)
    y(t)
    T(t)
  Dxt(t)
  Dyt(t)
 Dytt(t)
 Dxtt(t)

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

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

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

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

数値検索における数値的許容誤差を指定するオプション セットを作成します。

opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));

decic を使用して、変数とその導関数で整合性のある初期値を求めます。

[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt)
y0 =
    0.4771
   -0.8788
   -8.6214
         0
    0.0000
   -2.2333
   -4.1135

yp0 =
         0
    0.0000
         0
         0
   -2.2333
         0
         0

手順 6. ode15i を使用して DAE を求める

系を時間範囲 0t0.5 で積分して求解します。プロットにグリッド線と凡例を追加します。

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

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

legend(S, 'Location', 'Best')
grid on

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

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

r = 2;
F = @(t, Y, YP) f(t, Y, YP, g, m, r);

y0est = [r*cos(pi/6); r*sin(pi/6); 0; 0; 0; 0; 0];
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);

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

[tSol,y] = ode15i(F, [0, 0.5], y0, yp0, opt);
plot(tSol,y(:,1:origVars),'-o')

for k = 1:origVars
  S{k} = char(DAEvars(k));
end
legend(S, 'Location', 'Best')
grid on

関連するトピック