dsolve/numeric/DAE - 微分代数の初期値問題の数値解を探索
|
使い方
|
|
dsolve(daesys, numeric, vars, options)
|
|
パラメータ
|
|
daesys
|
-
|
セットまたはリスト、常微分方程式、代数方程式 (複数可)、および初期条件
|
numeric
|
-
|
数値解を検索するように dsolve に指示する名前
|
vars
|
-
|
(オプション) daesys に対する従属変数、あるいは従属変数のセットかリスト set or list of dependent variables
|
options
|
-
|
(オプション) keyword = value の形式の方程式
|
|
|
|
|
説明
|
|
•
|
numeric または type=numeric オプションを指定した dsolve コマンドおよび実数値による微分代数の初期値問題 (DAE IVP) は、DAE IVP の数値解を探索します。オプションの方程式 method=numericmethod が指定された場合 (この numericmethod は rkf45_dae、rosenbrock_dae または mebdfi のいずれか)、dsolve はこの方法を利用して数値解を取得します。
|
|
dsolve は、ほとんどの場合に問題が DAE 系であれば ODE 系とは対照的にこれを検出することができます。つまり、従属変数に対する純粋な代数方程式が存在する場合に検出が可能です。純粋な代数方程式を含まない DAE 系を代入しれた場合、これが DAE 系であることを指定するためにこの方法を 含めなければなりません。
|
•
|
制約された機械系で DAE の問題が発生することがよくあります (後述の振子の例を参照してください)。
|
'output'
|
=
|
キーワードまたは配列
|
'range'
|
=
|
数値..数値
|
'abserr'
|
=
|
数値またはリスト
|
'relerr'
|
=
|
数値
|
'stiff'
|
=
|
ブール式
|
'known'
|
=
|
名前または名前のリスト
|
'optimize'
|
=
|
ブール式
|
|
|
|
例外的に、DAE ソルバでは、コンポーネントリストごとに絶対的な誤差許容を指定することができます。予期される挙動については、問題の変数もリストとして指定しなければならず、'abserr' のエントリは、リストの変数に 1 対 1 で対応するか、系 (リストの変数の階数を使用して 1 階に変換される) の変数に 1 対 1 で対応する必要があります。詳細については、?dsolve[numeric,IVP] ヘルプページを参照してください。
|
•
|
デフォルトの DAE IVP 法は、改定 Runge-Kutta Fehlberg 法でなければなりません。これは 4-5 の基底階数の問題で使用する方法ですが、DAE 問題の解を探索できるように変更されています。デフォルトの stiff 方法は Rosenbrock 法で、この方法は、3-4 の基底階数を使用します。DAE 解を得るためのこれらの方法に加えられた拡張の詳細については、ヘルプページの ?dae_extension を参照してください。DAE IVP に使用可能な他の方法としては dsolve[mebdfi] 法があります。この名称は Modified Extended Backward-Differentiation Formula Implicit 方法の省略したものです。
|
•
|
一般に、DAE IVP ソルバは標準の微分 IVP ソルバに非常に似ています。そのため、このページでは主に両者の相違点について説明します。
|
•
|
現在、DAE ソルバは実数値問題の解の探索に用途が制限されています。
|
•
|
これらの方法を利用するには、系に隠れている制約事項を指定した初期条件ですべて満たす必要があります (つまり、DAE に関する初期条件に整合するセットでなければなりません)。これを満たしていない場合、計算結果がエラーになり、満たしていない条件についての情報が表示されます。
|
|
場合によっては、問題の初期条件が満たされるように fsolve で計算を行う必要があります。
|
•
|
?dae_extension 法の使用中は、微分のオプションは true に設定され、系は代数変数に対し十分に線形となるため (導関数のない変数が代入系に現れます)、これらの変数に対し、初期条件を省略することが可能です。必要な場合に初期条件を省略すると、エラーが生成されます。
|
•
|
以下のオプションも DAE 法で利用可能です。一部のオプションは DAE で利用できない場合があります。
|
'minstep'
|
=
|
数値
|
'maxstep'
|
=
|
数値
|
'initstep'
|
=
|
数値
|
'startinit'
|
=
|
ブール式
|
'events'
|
=
|
リスト
|
'event_pre'
|
=
|
キーワード
|
'event_maxiter'
|
=
|
整数
|
'projection'
|
=
|
ブール式
|
'differential'
|
=
|
ブール式
|
'implicit'
|
=
|
ブール式
|
'parameters'
|
=
|
名前のリスト
|
|
|
|
'minstep'、'maxstep'、および 'initstep'
|
|
デフォルトの方法である rkf45_dae および rosenbrock_dae では、'range' を指定すると、独立変数の新しい値が初期点と他の事前に要求されていた解点との間にない場合にのみ解が再計算されます。ここに積分の方向を一度も逆転しなかった効果が表れていて、すでに計算済みの間隔の解の評価もごく安価になります (解の値は補間で得られます)。解の保存も引数 storage を使用すると有効にできます。startinit パラメータでも、解の値が必要な場合はいつでもこの方法で再計算されるように強制できます。
|
|
'projection'、'differential'、および 'implicit'
|
|
展開法 に特有な 'projection'、'differential'、'implicit' の各オプションについて説明します。
|
•
|
指定の問題の代入値の解を計算することに加えて、プロシージャの戻り (output=procedurelist、listprocedure または operator) によって、parameters を指定する機能など、他の対話的な機能が追加されます。これらの機能の詳細については、dsolve[numeric,interactive] のページを参照してください。mebdfi 法はパラメータ化された問題を扱うことはできないことが例外です。
|
|
|
例
|
|
最初の例と同様に、ここでは 2-D デカルト座標の単位長さを持つ紐の質量の動力学をモデリングする問題 (振子の問題) を取り上げます。紐の質量位置を とし、速度を とします。
>
|
|
>
|
|
>
|
|
| (4.1) |
>
|
|
| (4.2) |
この質量は固定長 を持つ紐のため、制約は次のようになります。
>
|
|
| (4.3) |
ここで Euler-Lagrange 公式を使用して DAE 系の解を求めます。したがって、運動エネルギー とポテンシャルエネルギー を次のように計算します。
>
|
|
| (4.4) |
>
|
|
| (4.5) |
ここで は質量、 は動定数 (9.8 とします) です。Lagrangian 法と改定 Lagrangian 法で次のように指定します。
>
|
|
| (4.6) |
>
|
|
| (4.7) |
これで、Euler-Lagrange 公式が次のように構成できます。
>
|
|
| (4.8) |
>
|
|
| (4.9) |
>
|
|
| (4.10) |
>
|
|
| (4.11) |
これで振子の運動の方程式ができました。次に、整合する初期条件を決定する必要があります。このため、系の隠れた制約を見つけ出さなければなりません。制約は 1 つしかないため簡単に見つけることができます。
>
|
|
| (4.12) |
>
|
|
| (4.13) |
初期条件は、初期地点において 、、および を満たす必要があります。条件の自由度は 2 度だけです。したがって、最小値 で始動し、初期水平速度が の振子では、次のようになります。
>
|
|
| (4.14) |
>
|
|
| (4.15) |
これで上記のように、単位長さが で単位質量が 、で単位質量が の振子では、DAE 系と初期条件は次のようになると考えられます。
>
|
|
| (4.16) |
>
|
|
| (4.17) |
これで次のような解が得られます。
>
|
|
| (4.18) |
>
|
|
| (4.19) |
>
|
|
| (4.20) |
rosenbrock_dae による解は次のようになります。
>
|
|
| (4.21) |
>
|
|
| (4.22) |
>
|
|
| (4.23) |
mebdfi による解は次のようになります。
>
|
|
| (4.24) |
>
|
|
| (4.25) |
>
|
|
| (4.26) |
ここで、上記と同様の問題について考えますが、今度は、最初の紐の先に別の紐をつなげて第 2 の質量を追加します。追加する紐の長さは最初の紐の半分です (二重振子)。この系は次のように取得および求解できます。
>
|
|
| (4.27) |
>
|
|
| (4.28) |
>
|
|
| (4.29) |
第 2 の質量の曲線軌道は、次のようにしてプロットできます。
>
|
|
Warning, cannot evaluate the solution further right of 4.9606298, maxfun limit exceeded (see ?dsolve,maxfun for details)
| |
展開法を使用したため、 と に初期条件を指定していませんが、系はラムダに対して正しく線形となりました。
|
|
参照
|
|
dsolve/dae_extension, dsolve/Error_Control, dsolve[Events], dsolve[maxfun], dsolve[mebdfi], dsolve[numeric,interactive], dsolve[numeric,IVP], dsolve[numeric], dsolve[rkf45], dsolve[rosenbrock], dsolve[Stiffness], plots[odeplot]
|
|
Download Help Document
Was this information helpful?