dsolve/numeric/rosenbrock - 硬い常微分方程式の数値解を求める
使い方
dsolve(odesys, numeric, method=rosenbrock, vars, options)
dsolve(numeric, method=rosenbrock, procopts, options)
パラメータ
odesys - 集合またはリスト ; ( 連立 ) 常微分方程式および初期条件
numeric - 名前 ; dsolve に数値解を要求します
method=rosenbrock - method=rosenbrock を指定 ; 数値解法として使用します
vars - ( オプション ) 常微分方程式の未知変数を含む全ての一変数の不定関数, あるいはそれらの集合かリスト
options - ( オプション ) keyword = value 形式の方程式
procopts - 連立常微分方程式に用いる手続きの指定オプション (procedure, initial, start, number, procvars, jacobian のいずれか。 ) 。詳細は dsolve[numeric,IVP] をご参照ください。
|
説明
|
|
•
|
オプションに numeric および method=rosenbrock を指定した dsolve コマンドは、3 次の補間を含む 3-4 次 Runge-Kutta Rosenbrock 法 (Implicit Rosenbrock third-fourth order Runge-Kutta method) を用いて数値解を求めます。stiff を用いる場合、この解法は初期値問題を type=numeric と指定して解く場合のデフォルトの解法となっています。
|
|
Modes of Operation ( 操作方法? )
|
|
•
|
rosenbrock による解法には、二つの異なった ( 手続き型出力のための ) 操作方法があります。
|
|
range オプションを用いる場合、この解法は、指定範囲 (range) において、返された手続きを全て呼び出すために要求された解の値を高速に補間し、それを解の情報として内部に記憶して初期値問題 (IVP) の解を計算します。
|
|
この操作方法では、steppast が false としてデフォルト設定されており、積分の端点でステップが端点を越えてステップを踏まないように調整すべきであることを表示しています。指定範囲外の点で呼び出された手続きを返すことも可能ですが、あまり薦められません。
|
|
この方法は、適応型プロットとして odeplot の refine オプションと組み合わせて使用することができます ( つまり、refine が指定されたとき、odeplot はプロットのために前計算された点を使います。) 。
|
|
この方法は、解が特異になれるような問題に使用することは薦められません。各ステップを保存すると、特異性の近くでは多くのステップを取ることになり、メモリの使用量が深刻な問題となります。
|
|
この方法を使用している間に記憶する補間の情報は、以下で説明するようように interpolation=false オプションを用いることで無効にすることができます。これは、 ( 離散解に加えて ) 補間の記憶により多くのメモリを必要とする場所で高精度な解を得るために薦められます。解の値は最も近い 4 点の補間から得られるため、補間を無効にすることは一般的に薦められません。さらに、必ずしも 3 つのエラーのを補間に提供するわけではありません。
|
|
range オプションを用いない場合、初期値問題 (IVP) の解の値は記憶されませんが、要求される度に計算することになります。
|
|
この操作方法では、steppast が true としてデフォルト設定されており、解の値の出力は補間を使用する方法によって計算されます。このアプローチは、複数の点が必要となる場合により効率的な解の計算を提供しています。
|
|
この方法を含む補間の使用は、以下で説明する interpolation=false オプションを使用することによって無効にすることができますが、常に steppast=false を指定することが必要となります。
|
|
全ての解の値が記憶されるわけではないので、初期点と最も新しく計算された点の間の点が要求された場合はいつも、( 積分の反転を避けるために ) 初期値から計算しなおさなければなりません。したがって、初期値から移動した解の値を集めることが賢明です。
|
|
|
オプション
|
|
•
|
以下のオプションは rosenbrock で使用することが可能です。
|
'output' = keyword or array
'known' = name or list of names
'abserr' = numeric
'relerr' = numeric
'initstep' = numeric
'maxfun' = integer
'number' = integer
'procedure' = procedure
'jacobian' = procedure
'start' = numeric
'initial' = array
'procvars' = list
'startinit' = boolean
'implicit' = boolean
'optimize' = boolean
'range' = numeric..numeric
'interpolate' = boolean
'steppast' = boolean
'stop_cond' = list
'complex' = boolean
|
dsolve の出力形式を指定します。procedurelist, listprocedure, および operator は手続き型が出力を指定し、piecewise は指定された範囲の独立変数の値で piecewise 関数の形式で出力することを指定します。さらに 1-D 配列および配列は固定された独立変数の値を出力します。詳細は dsolve[numeric] をご参照ください。
|
|
ユーザが定義した既知関数を指定します。基本的な用法に関しては dsolve[numeric] をご参照ください。ここに書かれているように、numeric を指定した評価関数で定義された手続きに要件を加えます。もし numeric を指定しなかった場合は、`diff/` の規則を必要とします。
|
|
この `diff/` の規則は、独立変数、従変数および既知関数の計算に用いた導関数で表された関数の偏導関数 (partial derivatives) を示さなければなりません。注意 : ここでのキーワードは 偏微分 (partial) であり、その例は以下で見ることができます。これは taylorseries による計算法とは対照的なものです。独立変数による diff の規則は全微分 (total derivative) を示さなければなりません。
|
|
独立変数、従変数、および導関数で表された単独の偏導関数に関してだけは、各既知 (known) 関数を必要とします。これらは、`diff/` の規則を必要とするのではなく、それら自身の手続きとしてコード化することができます。
|
|
abserr, relerr, and initstep
|
|
解の精度と、解法の最初のステップ幅を指定します。詳細は dsolve[Error_Control] をご参照ください。rosenbrock でのデフォルトの値は互いに同じ値として abserr=Float(1,-6) と relerr=Float(1,-6) となっています。もし指定しなかった場合、initstep の値は連立常微分方程式の局所的な振る舞いを考慮して、解法が決定します。
|
|
一階の連立常微分方程式の右辺の評価数の最大値を指定します。このオプションは maxfun=0 と指定することで無効にできます。rosenbrock でのデフォルトの値は 30000 です。
|
|
number, procedure, start, initial, and procvars
|
|
このオプションは、procedure オプションおよびそれに関連するヤコビ行列 (Jacobian) の手続きの定義を指定するオプションと組み合わせて使用することができます。もし示されないならば、ヤコビ行列 (Jacobian) は入力した方程式系か、独立変数、従変数、あるいは結果を利用した入力手続きを使用して計算されます。ただし、このアプローチは入力の手続き形式の方程式が条件文、あるいはループを含んでいる場合は働かないことに注意してください。したがって、これらの場合では、手続き形式として jacobian を指定する必要があることがわかります。
|
|
手続き (procedure) を含むオプションとして、初期値問題の方程式系は最初に一階の方程式系として書かれなければなりません ( 詳細は dsolve[numeric,IVP] をご参照ください。 ) 。さらに、 jacobian の形式を得ることが出来ます。入力の手続きは、x を ( 入力した ) 独立変数の値、y を ( 入力した ) 従変数の値によるベクトル、fx を ( 出力した ) x を含む連立微分方程式の右辺の導関数の値のベクトル、および fy を ( 出力した ) y を含む各従変数連立常微分方程式の右辺の導関数の値の行列とした jac(x,y,fx,fy) の形式でなければなりません。
|
|
例として、いくつかの f(t,x(t),y(t)), g(t,x(t),y(t)) の形式を指定するために連立常微分方程式 diff(x(t),t) = f(t,x(t),y(t)), diff(y(t),t) = g(t,x(t),y(t)) を考えます。procedure と jacobian は以下のように手続きされます。 :
|
prc := proc(n,t,v,vp)
vp[1] := f(t,v[1],v[2]);
vp[2] := g(t,v[1],v[2]);
end proc:
jac := proc(t,v,ft,fv)
ft[1] := D[1](f)(t,v[1],v[2]); ft[2] := D[1](g)(t,v[1],v[2]);
fv[1,1] := D[2](f)(t,v[1],v[2]); fv[1,2] := D[3](f)(t,v[1],v[2]);
fv[2,1] := D[2](g)(t,v[1],v[2]); fv[2,2] := D[3](g)(t,v[1],v[2]);
end proc:
|
startinit,implicit, and optimize
|
|
解の値が必要とする独立変数の値の範囲を決定します。このオプションの使用は、dsolve[numeric] に説明がある手続き型の出力形式 (output types) のために計算法の振る舞いを大きく変えることになります (上記の range を使用するか否かについてをご参照ください。 ) 。
|
|
もし補間が積分のステップの間に加えられた解の値を計算すべきであるならば、論理演算子の値を受け入れます。これはデフォルトでは true と設定されています。ゆるい許容誤差 ( デフォルトと同様な ) のために、より計算された点、および滑らかなプロット、さらに補間された解のための一貫した精度を与えます。指定した範囲でより厳しい許容誤差を持つ問題には、この false 設定によって計算した解を記憶するメモリの量を減らすことを薦めます。
|
|
(range が指定されている場合に) 計算に範囲外の値を利用するか、(range が指定されていいない場合に) 要求した点で値を計算するかを計算法に指示した場合、論理演算子の値を受け入れます。デフォルトの値は range オプションが存在するかに依存します。範囲 (range) が指定されていれば false となり、それ以外では true となります。
|
|
記憶を行わない場合、このオプションは 補間 (interpolate) を含み、より一定の解の値を与え、解がより効果的に連続であるように同意しなければなりません。この場合、steppast=false, interpolate=false は計算法が出力点の近くでステップ幅を調整しなくてはならず、計算を続けるためにステップ幅を以前の調整に戻さなくてはならならことを意味します。
|
|
停止の評価基準のリストを受け入れます。このリストは、常微分方程式の解の数値積分をモニタする評価基準となります。評価基準の一つが満たされたとき、積分は停止し、現在のポイントを解として与えます。詳細は dsolve[Stop_Conditions] をご参照ください。
|
|
問題が複素数値である ( またはそうなるような ) 場合を表示する論理演算子を受け入れます。デフォルトでは方程式系と初期データの入力に基づいて検出するように設定されていますが、入力した方程式系が手続き型で定義されている場合や、初期データが実数である場合には、解を得るために complex=true を指定する必要があります。これは最初の実方程式系が複素方程式系になる可能性があります。この変化が起こる点は特異であると考えられるので、もし complex=true が指定されていなければ、この点で積分が停止することになります。
|
|
|
追加項目
|
|
•
|
infolevel[dsolve] を 2 と設定することで、最後の格子点での評価値の情報に誤差が含まれます。
|
•
|
rosenbrock による方法は、Digits の環境変数を変更して計算精度を上げることで、連立初期値問題 (IVPs) の解を高精度に計算することが可能になります。
|
|
|
|
例
|
|
-1, -100 の固有値を持つ硬い連立常微分方程式の区間 0..10 での強い解 :
>
|
deq1 := diff(y(t),t,t) + 101*diff(y(t),t) + 100*y(t)=0:
ic1 := y(0)=1, D(y)(0)=-1:
dsol1 := dsolve({deq1,ic1}, numeric, stiff=true,
range=0..10);
|
| (2.1) |
| (2.2) |
原子炉動特性に見られる硬い非線形方程式系の解の値の計算 :
>
|
deq2 := diff(u(t),t) = 0.01 - (1+(u(t)+1000)*(u(t)+1))*(0.01+u(t)+v(t)),
diff(v(t),t) = 0.01 - (1+v(t)^2)*(0.01+u(t)+v(t));
|
| (2.3) |
>
|
ic2 := u(0)=0, v(0)=0:
dsol2 := dsolve({deq2,ic2}, numeric, method=rosenbrock);
|
| (2.4) |
| (2.5) |
| (2.6) |
( diff の規則を必要とした ) 多重関数に対して known オプションを使用した例。以下の関数をしようします。
>
|
f(x,y(x)) = Int(exp(-t*2/100)/10,t=0..x)-y(x);
|
| (2.7) |
'f' の定義
>
|
f := proc(x,y) local t;
if not type(evalf(x),numeric) or
not type(evalf(y),numeric) then
'procname'(x,y);
else
evalf(Int(exp(-t*2/100)/10,t=0..x))-y;
end if;
end proc;
|
| (2.8) |
さらに x と y を含む偏導関数に `diff/` の規則を使用します。 :
>
|
`diff/f` := proc(x,y);
if args[3]=x then
diff(x,args[3])*exp(-x*2/100)/10;
elif args[3]=y then
-1;
else
error "unable to differentiate";
end if;
end proc;
|
| (2.9) |
この情報を用いれば、解の手続きを得るために dsolve を呼び出すことで、すぐに解の値を得ることができます。
>
|
dsol := dsolve({diff(y(x),x)=f(x,y(x)), y(0)=0}, numeric,
known=f,method=rosenbrock);
|
| (2.10) |
| (2.11) |
|
|
参照
|
|
dsolve[classical], dsolve[dverk78], dsolve[Error_Control], dsolve[gear], dsolve[lsode], dsolve[maxfun], dsolve[numeric], dsolve[numeric,IVP], dsolve[rkf45], dsolve[Stiffness], dsolve[Stop_Conditions], dsolve[taylorseries], infolevel, plots[odeplot]
|
|
参考文献
|
|
|
Hairer, E. and Wanner, G. "Solving Ordinary Differential Equations II, 2nd ed." Springer (1996).
|
|
Shampine, L.F. and Corless, R.M. "Initial Value Problems for ODEs in Problem Solving Environments." J. Comp. Appl. Math, vol. 125(1-2) (2000) 31-40.
|
|
|