24.2 Differential-Algebraic Equations

Last-modified: 2025-03-07 (金) 22:56:39

24.2 微分代数方程式

この関数は、daspk次の形式の微分代数方程式を解くために使用できる。

0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0

ここで、 x-dotはx の導関数です。この方程式は、 Petzold の DAE ソルバーDASPKを使用して解きます。

: daspk [x, xdot, istate, msg] = (fcn, x_0, xdot_0, t, t_crit)

一連の微分代数方程式を解きます。

daspk方程式を解く

0 = f (x, xdot, t)

x(t_0) = x_0、xdot(t_0) = xdot_0

解は行列xとxdotで返され、結果行列の各行はベクトルtの要素の 1 つに対応します。 tの最初の要素 はt_0で、システムの初期状態x_0とその導関数xdot_0に対応します。そのため、出力xの最初の行はx_0で、出力xdotの最初の行はxdot_0になります。

最初の引数fcnは、方程式の集合の残差ベクトルを計算するために呼び出す 関数fを指定する文字列、インライン、または関数ハンドルです。次の形式である必要があります。

res = f ( x、xdot、t )

です。
ここで、 x、xdot、res はベクトルであり、tはスカラーです。

fcnが2要素の文字列配列または文字列、インライン関数、または関数ハンドルの2要素のセル配列である場合、最初の要素は上記の関数fの名前であり、2番目の要素は修正ヤコビアンを計算する関数の名前です。

      df       df
jac = -- + c ------
      dx     d xdot

修正ヤコビ関数は次の形式を持つ必要がある。

jac = j (x, xdot, t, c)

2 番目と 3 番目の引数はdaspk状態とその導関数の初期条件を指定し、4 番目の引数は初期条件に対応する時間を含む、ソリューションが求められる出力時間のベクトルを指定します。

初期状態と導関数のセットは、必ずしも一貫している必要はありません。 一貫していない場合は、関数を使用して 追加情報を提供し、一貫した開始点を計算できる daspk_optionsようにする必要があります。daspk

5 番目の引数はオプションであり、DAE ソルバーが積分しない過去の時間のセットを指定するために使用できます。これは、特異点や導関数に不連続がある点に関する問題を回避するのに役立ちます。

計算が成功すると、istateの値はゼロより大きくなります ( DASPKの Fortran バージョンと一致します)。

計算が成功しなかった場合、istateの値はゼロ未満になり、msg には追加情報が含まれます。

関数を使用してdaspk_options、 のオプション パラメータを設定できますdaspk。

See also: dassl.

: daspk_options ()

: val = daspk_options (opt)

: daspk_options (opt, val)

daspk関数のオプションを照会または設定します。

引数なしで呼び出されると、使用可能なすべてのオプションの名前とその現在の値が表示されます。

引数を 1 つ指定すると、オプションoptの値を返します。

2 つの引数で呼び出されると、daspk_optionsオプション opt が値valに設定されます。

オプションには以下が含まれます

"absolute tolerance"

絶対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、相対許容差も同じ長さのベクトルである必要があります。

"relative tolerance"

相対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、絶対許容差も同じ長さのベクトルである必要があります。

各積分ステップで適用されるローカルエラーテストは

abs (local error in x(i))
      <= rtol(i) * abs (Y(i)) + atol(i)

"compute consistent initial condition"
状態ベクトルの微分変数を「Y_d' および代数変数は 'やあ' は、ddaspk次の 2 つの初期化問題のいずれかを解決できます。

Y_dが与えられたら、Y_aとY'_dを計算する
Y'が与えられたら、Yを計算します。

どちらの場合も、指定されたコンポーネントの初期値が入力され、未知のコンポーネントの初期推定値も入力として提供する必要があります。最初の問題を解決するにはこのオプションを 1 に設定し、2 番目の問題を解決するには 2 に設定します (デフォルトは 0 であるため、一貫した初期条件のセットを提供する必要があります)。

このオプションをゼロ以外の値に設定する場合は、 "algebraic variables"問題内のどの変数が代数的であるかを宣言するオプションも設定する必要があります。

"use initial condition heuristics"

以下で説明する初期条件ヒューリスティック オプションを使用するには、ゼロ以外の値に設定します。

"initial condition heuristics"

初期条件の計算を制御するために使用できる次のパラメータのベクトル。

MXNIT

ニュートン反復の最大回数 (デフォルトは 5)。

MXNJ

ヤコビ評価の最大数 (デフォルトは 6)。

MXNH

オプションが 1 に設定されている場合、試行される人工的なステップサイズ パラメータの値の最大数"compute consistent initial condition"(デフォルトは 5)。

許可されるニュートン反復の最大合計数は、 オプションが 1 に設定されているMXNIT*MXNJ*MXNH場合と2 に設定されている場合であることに注意してください。 "compute consistent initial condition"MXNIT*MXNJ

LSOFF

ラインサーチアルゴリズムを無効にするには、ゼロ以外の値に設定します (デフォルトは 0)。

STPTOL

ラインサーチアルゴリズムにおける最小スケールステップ(デフォルトはeps^(2/3))。

EPINIT

ニュートン反復収束テストにおけるスイング係数。テストは、近似ヤコビアンで事前に乗算された残差ベクトルに適用されます。収束するには、このベクトルの重み付き RMS ノルム (誤差の重みでスケール) が 未満である必要がありますEPINIT*EPCON。ここで、EPCON= 0.33 は、時間ステップで使用される類似のテスト定数です。デフォルトはEPINIT= 0.01 です。

"print initial condition info"

初期条件計算に関する詳細情報を表示するには、このオプションを 0 以外の値に設定します (デフォルトは 0)。

"exclude algebraic variables from error test"

代数変数をエラー テストから除外するには、ゼロ以外の値に設定します。"algebraic variables"問題内のどの変数が代数変数であるかを宣言するオプションも設定する必要があります (デフォルトは 0)。

"algebraic variables"

状態ベクトルと同じ長さのベクトル。非ゼロ要素は、状態ベクトルの対応する要素が代数変数であることを示します (つまり、その導関数は方程式セットに明示的に表示されません)。

"compute consistent initial condition"このオプションは、および オプションで必須です "exclude algebraic variables from error test"。

"enforce inequality constraints"

オプションで指定された不等式制約を適用するには、次のいずれかの値に設定します"inequality constraint types" (デフォルトは 0)。

初期条件の計算でのみ制約チェックを実行します。
統合中に制約チェックを強制します。
オプション 1 と 2 の両方を適用します。
"inequality constraint types"
不等式制約のタイプを指定する状態と同じ長さのベクトル。ベクトルの各要素は状態の要素に対応し、次のコードのいずれかを割り当てる必要があります。

-2
ゼロ未満。
-1
ゼロ以下。
0
制約はありません。
1
ゼロ以上。
2
ゼロより大きい。

このオプションは、オプションがゼロ以外の場合にのみ効果があります "enforce inequality constraints"。

"initial step size"

微分代数問題では、最初のステップで深刻なスケーリングの問題が発生する場合があります。問題のスケーリングについて十分な知識がある場合は、初期ステップ サイズ (既定値は自動的に計算されます) を指定して、この問題を軽減することができます。

"maximum order"

解法の最大次数を制限します。このオプションは 1 から 5 までの範囲で指定する必要があります (デフォルトは 5)。

"maximum step size"

最大ステップサイズを設定すると、非常に大きな領域を通過することを回避できます (デフォルトは指定されていません)。

Octave には、制約 (停止条件) 付きの DAE を解くために使用できる 、 DASPKの以前のバージョンであるDASSLとDASRTも含まれています。

: [x, xdot, istate, msg] = dassl (fcn, x_0, xdot_0, t, t_crit)

一連の微分代数方程式を解きます。

dassl方程式を解く

0 = f (x, xdot, t)

x(t_0) = x_0、xdot(t_0) = xdot_0

解は行列xとxdotで返され、結果行列の各行はベクトルtの要素の 1 つに対応します。 tの最初の要素 はt_0で、システムの初期状態x_0とその導関数xdot_0に対応します。そのため、出力xの最初の行はx_0で、出力xdotの最初の行はxdot_0になります。

最初の引数fcnは、方程式の集合の残差ベクトルを計算するために呼び出す 関数fを指定する文字列、インライン、または関数ハンドルです。次の形式である必要があります。

res = f (x, xdot, t)

です。
ここで、 x、xdot、res はベクトルであり、tはスカラーです。

fcnが2要素の文字列配列または文字列、インライン関数、または関数ハンドルの2要素のセル配列である場合、最初の要素は上記の関数fの名前であり、2番目の要素は修正ヤコビアンを計算する関数の名前です。

      df       df
jac = -- + c ------
      dx     d xdot

修正ヤコビ関数は次の形式を持つ必要がある。

jac = j (x, xdot, t, c)

2 番目と 3 番目の引数はdassl状態とその導関数の初期条件を指定し、4 番目の引数は初期条件に対応する時間を含む、ソリューションが求められる出力時間のベクトルを指定します。

初期状態と導関数のセットは、厳密には一貫している必要はありません。ただし、実際には、DASSL は一貫性のあるセットを決定するのがあまり得意ではないため、初期値によって関数がゼロに評価されるようにすることが最善です。

5 番目の引数はオプションであり、DAE ソルバーが積分しない過去の時間のセットを指定するために使用できます。これは、特異点や導関数に不連続がある点に関する問題を回避するのに役立ちます。

計算が成功すると、istateの値はゼロより大きくなります ( DASSLの Fortran バージョンと一致します)。

計算が成功しなかった場合、istateの値はゼロ未満になり、msg には追加情報が含まれます。

関数を使用してdassl_options、 のオプション パラメータを設定できますdassl。

See also: daspk, dasrt, lsode.

: dassl_options ()

: val = dassl_options (opt)

: dassl_options (opt, val)

dassl関数のオプションを照会または設定します。

引数なしで呼び出されると、使用可能なすべてのオプションの名前とその現在の値が表示されます。

引数を 1 つ指定すると、オプションoptの値を返します。

2 つの引数で呼び出されると、dassl_optionsオプション opt が値valに設定されます。

オプションには以下が含まれます

"absolute tolerance"

絶対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、相対許容差も同じ長さのベクトルである必要があります。

"relative tolerance"

相対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、絶対許容差も同じ長さのベクトルである必要があります。

各積分ステップで適用されるローカルエラーテストは

 abs (local error in x(i)) <= ...
     rtol(i) * abs (Y(i)) + atol(i)
"compute consistent initial condition"

ゼロ以外の場合、dassl一貫した初期条件セットを計算しようとします。これは通常信頼できないため、一貫したセットを提供し、このオプションをゼロに設定しておくのが最善です。

"enforce nonnegativity constraints"

方程式の解が常に非負になることが分かっている場合は、このパラメータをゼロ以外の値に設定すると役立つ場合があります。ただし、最初はこのオプションをゼロに設定したままにしておき、うまく機能しない場合にのみゼロ以外の値に設定するのがおそらく最善です。

"initial step size"

微分代数問題では、最初のステップで深刻なスケーリングの問題が発生する場合があります。問題のスケーリングについて十分な知識がある場合は、初期ステップ サイズを指定してこの問題を軽減することができます。

"maximum order"

解法の最大次数を制限します。このオプションは 1 から 5 までの範囲でなければなりません。

"maximum step size"

最大ステップサイズを設定すると、非常に大きな領域を通過することを回避できます (デフォルトは指定されていません)。

"step limit"

基礎となる Fortran コードへの 1 回の呼び出しで試行する統合ステップの最大数。

: [x, xdot, t_out, istat, msg] = dasrt (fcn, g, x_0, xdot_0, t)

: … = dasrt (fcn, g, x_0, xdot_0, t, t_crit)

: … = dasrt (fcn, x_0, xdot_0, t)

: … = dasrt (fcn, x_0, xdot_0, t, t_crit)

一連の微分代数方程式を解きます。

dasrt方程式を解く

0 = f (x, xdot, t)

x(t_0) = x_0、xdot(t_0) = xdot_0

機能的な停止基準(ルート解決)付き。

解は行列xとxdotで返され、結果行列の各行はベクトルt_outの要素の 1 つに対応します。 tの最初の要素 はt_0で、システムの初期状態x_0とその導関数xdot_0に対応します。そのため、出力xの最初の行はx_0で、出力xdotの最初の行はxdot_0になります。

ベクトルt は積分の長さの上限を提供します。停止条件が満たされると、ベクトル t_out はtより短くなり、 t_outの最後の要素は停止条件が満たされたポイントになり、ベクトルtのどの要素とも一致しない場合があります。

最初の引数fcnは、方程式の集合の残差ベクトルを計算するために呼び出す 関数fを指定する文字列、インライン、または関数ハンドルです。次の形式である必要があります。

res = f (x, xdot, t)

です。
ここで、 x、xdot、res はベクトルであり、tはスカラーです。

fcnが2要素の文字列配列または文字列、インライン関数、または関数ハンドルの2要素のセル配列である場合、最初の要素は上記の関数fの名前であり、2番目の要素は修正ヤコビアンを計算する関数の名前です。

      df       df
jac = -- + c ------
      dx     d xdot

修正ヤコビ関数は次の形式を持つ必要がある。

jac = j (x, xdot, t, c)

オプションの2番目の引数は、積分中に根が求められる制約関数を定義する関数の名前です。この関数は次の形式である必要があります。

g_out = g ( x , t )

制約関数の値のベクトルを返します。制約関数のいずれかの値の符号が変わると、DASRT は 符号が変わった時点で積分を停止しようとします。

制約関数の名前を省略すると、daspkdasslまたはdasrtと同じ問題を解きます。

制約関数では丸めや積分誤差による数値誤差が発生するため、DASRT は誤った根を返す場合や、 Tのほぼ等しい値が 2 つ以上ある場合に同じ根を返す場合がある ことに注意してください。このような誤った根が疑われる場合は、制約関数の評価において、誤差許容値を小さくするか、精度を高くすることを検討する必要があります。

ある制約関数の根が問題の終わりを定義する場合、DASRTへの入力では、その根をわずかに超えた点までの積分を許可する必要があります。これにより、DASRT は補間によって根を特定できます。

3 番目と 4 番目の引数はdasrt、状態とその導関数の初期条件を指定します。4 番目の引数は、初期条件に対応する時間を含む、ソリューションが求められる出力時間のベクトルを指定します。

初期状態と導関数のセットは、厳密には一貫している必要はありません。ただし、実際には、DASSL は一貫性のあるセットを決定するのがあまり得意ではないため、初期値によって関数がゼロに評価されるようにすることが最善です。

6 番目の引数はオプションであり、DAE ソルバーが積分しない過去の時間のセットを指定するために使用できます。これは、特異点や導関数に不連続がある点に関する問題を回避するのに役立ちます。

計算が成功すると、istateの値はゼロより大きくなります ( DASSLの Fortran バージョンと一致します)。

計算が成功しなかった場合、istateの値はゼロ未満になり、msg には追加情報が含まれます。

dasrt関数を使用してdasrt_options、 のオプション パラメータを設定できます。

参照: dasrt_options、daspk、dasrt、lsode。

: dasrt_options ()
: val = dasrt_options (opt)
: dasrt_options (opt, val)
関数のオプションを照会または設定しますdasrt。

引数なしで呼び出されると、使用可能なすべてのオプションの名前とその現在の値が表示されます。

引数を 1 つ指定すると、オプションoptの値を返します。

2 つの引数で呼び出されると、dasrt_optionsオプション opt が値valに設定されます。

オプションには以下が含まれます

"absolute tolerance"

絶対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、相対許容差も同じ長さのベクトルである必要があります。

"relative tolerance"

相対許容差。ベクトルまたはスカラーのいずれかになります。ベクトルの場合は、状態ベクトルの次元と一致する必要があり、絶対許容差も同じ長さのベクトルである必要があります。

各積分ステップで適用されるローカルエラーテストは

 abs (x(i) のローカル誤差) <= ...
     rtol(i) * abs (Y(i)) + atol(i)
"initial step size"

微分代数問題では、最初のステップで深刻なスケーリングの問題が発生する場合があります。問題のスケーリングについて十分な知識がある場合は、初期ステップ サイズを指定してこの問題を軽減することができます。

"maximum order"

解法の最大次数を制限します。このオプションは 1 から 5 までの範囲でなければなりません。

"maximum step size"

最大ステップサイズを設定すると、非常に大きな領域を通過することを回避できます。

"step limit"

基礎となる Fortran コードへの 1 回の呼び出しで試行する統合ステップの最大数。

See K. E. Brenan, et al., Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, North-Holland (1989), DOI: https://doi.org/10.1137/1.9781611971224, for more information about the implementation of DASSL.