Main Content

decic

代数制約のある 1 階の陰的 ODE 系に整合する初期条件を求める

ページ内をすべて折りたたむ

構文

[y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)

説明

例

[y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options) は、関数 reduceDAEToODE によって返される代数制約を使用して陰的な 1 階常微分方程式系に整合する初期条件を求めます。

[eqs,constraintEqs] = reduceDAEToODE(DA_eqs,vars) の呼び出しにより、微分代数方程式系 DA_eqs は陰的 ODE 系 eqs に簡約されます。また、系の簡約化時に発生した拘束方程式を返します。この ODE 系の変数およびその導関数について、decic は時間 t0 における整合する初期条件 y0 と yp0 を求めます。

数値 y0、yp0 を微分方程式 subs(eqs, [t; vars(t); diff(vars(t))], [t0; y0; yp0]) および拘束方程式 subs(constr, [t; vars(t); diff(vars(t))], [t0; y0; yp0]) に代入するとゼロ ベクトルが生成されます。ここで、vars は列ベクトルでなければなりません。

y0_est は時間 t0 における変数 vars の数値推定を指定し、fixedVars は数値検索中に変化してはならない y0_est の値を示します。オプションの引数 yp0_est を使用すると、時間 t0 における変数 vars の導関数の値への数値近似を指定できます。

例

ODE 系に整合する初期条件を求める

DAE 系を陰的な ODE 系に簡約します。次に、得られた ODE 系の変数およびそれらの 1 次導関数に整合する初期条件を求めます。

次の微分代数方程式系を作成します。

syms x(t) y(t)
DA_eqs = [diff(x(t),t) == cos(t) + y(t),...
          x(t)^2 + y(t)^2 == 1];
vars = [x(t); y(t)];

reduceDAEToODE を使用してこの系を陰的 ODE 系に変換します。

[eqs, constraintEqs] = reduceDAEToODE(DA_eqs, vars)
eqs =
                 diff(x(t), t) - y(t) - cos(t)
 - 2*x(t)*diff(x(t), t) - 2*y(t)*diff(y(t), t)
 
constraintEqs =
1 - y(t)^2 - x(t)^2

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

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

時間の値 t0 = 0 と、変数およびその導関数の整合する値の数値推定を決定します。

t0 = 0;
y0_est = [0.1, 0.9];
yp0_est = [0.0, 0.0];

制約を、固定パラメーター y をもつ変数 x の代数方程式として扱うことができます。ここでは、fixedVars = [0 1] と設定します。または、固定パラメーター x をもつ変数 y の代数方程式として扱うこともできます。ここでは、fixedVars = [1 0] と設定します。

まず、初期値 x(t0) = y0_est(1) = 0.1 を設定します。

fixedVars = [1 0];
[y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)
y0 =
    0.1000
    0.9950

yp0 =
    1.9950
   -0.2005

ここで、fixedVars を [0 1] に変更します。これにより y(t0) = y0_est(2) = 0.9 が決まります。

fixedVars = [0 1];
[y0,yp0] = decic(eqs,vars,constraintEqs,t0,y0_est,fixedVars,yp0_est,options)
y0 =
   -0.4359
    0.9000

yp0 =
    1.9000
    0.9202

これらの初期値が方程式および制約を満たす整合性のある初期値であることを検証します。

subs(eqs, [t; vars; diff(vars,t)], [t0; y0; yp0])
ans =
 0
 0
subs(constraintEqs, [t; vars; diff(vars,t)], [t0; y0; yp0])
ans =
0

入力引数

すべて折りたたむ

陰的な常微分方程式系。シンボリック方程式またはシンボリック式のベクトルとして指定します。この式は右辺がゼロの方程式を表します。

通常は、reduceDAEToODE によって返された式を使用します。

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

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

系の簡約化時に発生する拘束方程式。シンボリック方程式または式のベクトルとして指定します。これらの式または方程式は変数 vars に依存しますが、その導関数には依存しません。

通常は、reduceDAEToODE によって返された拘束方程式を使用します。

初期時間。数値として指定します。

初期時間 t0 における変数 vars の値の推定。数値ベクトルとして指定します。

y0_est のどの要素が固定値であるかを示す入力ベクトル。0 または 1 をもつベクトルとして指定します。y0_est の固定値は fixedVars の値 1 に対応します。これらの値は数値検索中には変更されません。fixedVars のゼロは、y0_est の変数内の、decic が拘束方程式を解いて求める対象に対応します。0 の数は拘束方程式の数と一致しなければなりません。変数 vars(fixedVars == 0) に対し、この拘束方程式のヤコビ行列は可逆でなければなりません。

初期時間 t0 における変数 vars の 1 次導関数の値の推定。数値ベクトルとして指定します。

数値検索のオプション。odeset によって返されるオプション構造体として指定します。たとえば、ここで数値検索の許容誤差を指定できます。

出力引数

すべて折りたたむ

各変数について整合する初期値。数値列ベクトルとして返されます。

数値列ベクトルとして返される、各変数の 1 次導関数に対応する初期値。

バージョン履歴

R2014b で導入

参考

| | | | | | | | |

トピック