微分代数の初期値問題の数値解を探索

はじめる前に
新機能一覧
Maple ワークシートの作成
Mapleワークシートを共有
Maple ウィンドウのカスタマイズ
Maple システムのカスタマイズ
コネクティビティ
Mathematics
数学
- テンソル解析
- パッケージ
- ファンクションアドバイザ
- ベキ級数
- ベクトル解析
- 不活性関数
- 代数
- 因数分解と方程式解法
  - LinearFunctionalSystems
  - LinearOperators
  - LREtools
  - Numerical Solutions
    - fsolve
    - dsolve/DAE
    - dsolve/dverk78
    - dsolve/numeric
    - RootFinding
    - solve/float
  - QDifferenceEquations
  - RegularChains
  - SolveTools
  - 根
  - 解く
  - 記号解
  - hermite pade コマンド
  - invfunc
  - isolate
  - LeastSquares
  - LinearSolve
  - MinimalPolynomial コマンド
  - minimax コマンド
  - msolve
  - powsolve
  - rsolve
  - singular
- 基本数学
- 基本的な注意事項
- 変分法
- 変換
- 幾何
- 微分代数方程式
  - DAE 例題ワークシート
  - mebdfi オプション
  - 微分代数方程式（DAE）
  - 微分代数方程式の拡張
- 微分幾何学
- 微分方程式
- 微積分
- 数
- 数値計算
- 数学関数
- 数論
- 最適化
- 特殊関数
- 統計
- 線形代数
- 群論
- 評価
- 論理
- 過去のパッケージ
- 金融
- 離散数学
- piecewise examples
- グラフ電卓
- 区分関数
物理
プログラミング
グラフィックス
Student パッケージ
Science and Engineering
科学・エンジニアリング
アプリケーションと例題
リファレンス
システム
MapleSim
MapleSim ツールボックス
MapleSim ヘルプ
Tasks
Toolboxes
エラーメッセージガイド
マニュアル
数学アプリ

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 の問題が発生することがよくあります (後述の振子の例を参照してください)。

•	dsolve の戻り値、および以下の高レベルオプションや一般的なオプションの詳細については、dsolve[numeric] および dsolve[Error_Control] を参照してください。

'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'

これらのオプションは、dsolve[mebdfi] 法でのみ利用できます。詳細については、dsolve[Error_Control] を参照してください。

'startinit'= boolean

このオプションは、後続する計算に関する数値積分の挙動を制御します。詳細については、dsolve[numeric,IVP] を参照してください。

デフォルトの方法である rkf45_dae および rosenbrock_dae では、'range' を指定すると、独立変数の新しい値が初期点と他の事前に要求されていた解点との間にない場合にのみ解が再計算されます。ここに積分の方向を一度も逆転しなかった効果が表れていて、すでに計算済みの間隔の解の評価もごく安価になります (解の値は補間で得られます)。解の保存も引数 storage を使用すると有効にできます。startinit パラメータでも、解の値が必要な場合はいつでもこの方法で再計算されるように強制できます。

'events'= list

'event_pre'= keyword

'event_maxiter'= integer

これらのオプションは rkf45_dae 法と rosenbrock_dae 法でのみ利用可能で、イベントの挙動を指定および制御するために使用します。これらのオプションは通常の IVP 問題のオプションと同様です。詳細については、ヘルプページの ?dsolve[numeric,IVP] および ?dsolve[Events] を参照してください。

'projection'、'differential'、および 'implicit'

展開法に特有な 'projection'、'differential'、'implicit' の各オプションについて説明します。

•	指定の問題の代入値の解を計算することに加えて、プロシージャの戻り (output=procedurelist、listprocedure または operator) によって、parameters を指定する機能など、他の対話的な機能が追加されます。これらの機能の詳細については、dsolve[numeric,interactive] のページを参照してください。mebdfi 法はパラメータ化された問題を扱うことはできないことが例外です。

例

最初の例と同様に、ここでは 2-D デカルト座標の単位長さを持つ紐の質量の動力学をモデリングする問題 (振子の問題) を取り上げます。紐の質量位置をとし、速度をとします。

r := Vector(2, {(1) = x(t), (2) = y(t)})

(4.1)

v := Vector(2, {(1) = diff(x(t), t), (2) = diff(y(t), t)})

(4.2)

この質量は固定長を持つ紐のため、制約は次のようになります。

(4.3)

ここで Euler-Lagrange 公式を使用して DAE 系の解を求めます。したがって、運動エネルギーとポテンシャルエネルギーを次のように計算します。

(4.4)

(4.5)

ここでは質量、は動定数 (9.8 とします) です。Lagrangian 法と改定 Lagrangian 法で次のように指定します。

(4.6)

(4.7)

これで、Euler-Lagrange 公式が次のように構成できます。

$Q := {-2*lambda(t)*x(t)-m*(diff(x(t), `$`(t, 2))), -x(t)^2-y(t)^2+l^2, -m*g-2*lambda(t)*y(t)-m*(diff(y(t), `$`(t, 2)))}$

(4.8)

EQx := isolate(select(has, Q, diff(x(t), t, t))[1], diff(x(t), t, t))

EQx := diff(x(t), `$`(t, 2)) = -2*lambda(t)*x(t)/m

(4.9)

EQy := isolate(select(has, Q, diff(y(t), t, t))[1], diff(y(t), t, t))

EQy := diff(y(t), `$`(t, 2)) = -(m*g+2*lambda(t)*y(t))/m

(4.10)

(4.11)

これで振子の運動の方程式ができました。次に、整合する初期条件を決定する必要があります。このため、系の隠れた制約を見つけ出さなければなりません。制約は 1 つしかないため簡単に見つけることができます。

(4.12)

(4.13)

初期条件は、初期地点において、、およびを満たす必要があります。条件の自由度は 2 度だけです。したがって、最小値で始動し、初期水平速度がの振子では、次のようになります。

$sys := {-2*x(0)*(D(x))(0)-2*y(0)*(D(y))(0), -x(0)^2-y(0)^2+l^2, -2*(D(x))(0)^2+4*x(0)^2*lambda(0)/m-2*(D(y))(0)^2+2*y(0)*(m*g+2*lambda(0)*y(0))/m, y(0) = -l, (D(x))(0) = vx}$

(4.14)

$ini := {lambda(0) = (1/2)*m*(vx^2+l*g)/l^2, x(0) = 0, y(0) = -l, (D(x))(0) = vx, (D(y))(0) = 0}$

(4.15)

これで上記のように、単位長さがで単位質量が、で単位質量がの振子では、DAE 系と初期条件は次のようになると考えられます。

dsys := eval({EQc, EQx, EQy}, [l = 1, m = 1, g = 9.8, vx = 1/10])

$dsys := {-x(t)^2-y(t)^2+1, diff(x(t), `$`(t, 2)) = -2*lambda(t)*x(t), diff(y(t), `$`(t, 2)) = -9.8-2*lambda(t)*y(t)}$

(4.16)

dini := eval(ini, [l = 1, m = 1, g = 9.8, vx = 1/10])

(4.17)

これで次のような解が得られます。

(4.18)

[t = .500000000000000, lambda(t) = 4.89750023744205887, x(t) = 0.319392547405407887e-1, diff(x(t), t) = 0.564531449666508191e-3, y(t) = -.999489811893806701, diff(y(t), t) = 0.180382736946334369e-4]

(4.19)

[t = 1., lambda(t) = 4.90499904949563259, x(t) = 0.360893004470803490e-3, diff(x(t), t) = -0.999936779200488873e-1, y(t) = -.99999993487231431, diff(y(t), t) = -0.360870084253994927e-4]

(4.20)

rosenbrock_dae による解は次のようになります。

(4.21)

[t = .500000000000000, lambda(t) = 4.89750019118586088, x(t) = 0.319392364763693773e-1, diff(x(t), t) = 0.564626004686089771e-3, y(t) = -.999489813543333439, diff(y(t), t) = 0.180482728717063205e-4]

(4.22)

[t = 1., lambda(t) = 4.90499906207810721, x(t) = 0.360968463263104937e-3, diff(x(t), t) = -0.999936779452163804e-1, y(t) = -.99999993669131815, diff(y(t), t) = -0.360897908606885759e-4]

(4.23)

mebdfi による解は次のようになります。

(4.24)

[t = .500000000000000, lambda(t) = 4.89750245194350509, x(t) = 0.319388403694855147e-1, diff(x(t), t) = 0.564871203422257005e-3, y(t) = -.999489825189509129, diff(y(t), t) = 0.180670454630691393e-4]

(4.25)

[t = 1., lambda(t) = 4.90497944724429935, x(t) = 0.361097312194549836e-3, diff(x(t), t) = -0.999925039194764775e-1, y(t) = -.99999993426187950, diff(y(t), t) = -0.359881610159317974e-4]

(4.26)

ここで、上記と同様の問題について考えますが、今度は、最初の紐の先に別の紐をつなげて第 2 の質量を追加します。追加する紐の長さは最初の紐の半分です (二重振子)。この系は次のように取得および求解できます。

>	$dsysd := {diff(x2(t), t, t)-2lambda2(t)(x1(t)-x2(t)), (x1(t)-x2(t))^2+(y1(t)-y2(t))^2-1/4, x1(t)^2+y1(t)^2-1, diff(x1(t), t, t)+2lambda1(t)x1(t)+2lambda2(t)(x1(t)-x2(t)), diff(y2(t), t, t)+9.8-2lambda2(t)(y1(t)-y2(t)), diff(y1(t), t, t)+9.8+2lambda1(t)y1(t)+2lambda2(t)(y1(t)-y2(t))}$

$dsysd := {diff(x2(t), `$`(t, 2))-2*lambda2(t)*(x1(t)-x2(t)), (x1(t)-x2(t))^2+(y1(t)-y2(t))^2-1/4, x1(t)^2+y1(t)^2-1, diff(x1(t), `$`(t, 2))+2*lambda1(t)*x1(t)+2*lambda2(t)*(x1(t)-x2(t)), diff(y2(t), `$`(t, 2))+9.8-2*lambda2(t)*(y1(t)-y2(t)), diff(y1(t), `$`(t, 2))+9.8+2*lambda1(t)*y1(t)+2*lambda2(t)*(y1(t)-y2(t))}$

(4.27)

ics := {x1(0) = 0, x2(0) = 0, y1(0) = -1, y2(0) = -3/2, (D(x1))(0) = -3, (D(x2))(0) = 4, (D(y1))(0) = 0, (D(y2))(0) = 0}

(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]