dsolve/numeric - 常微分方程式の数値解を求める
|
使い方
|
|
dsolve(odesys, numeric, vars, options)
dsolve(numeric, procopts, options)
|
|
パラメータ
|
|
odesys
|
-
|
集合またはリスト; (連立) 常微分または微分代数方程式および初期/境界条件
|
numeric
|
-
|
名前; dsolve に数値解を求める命令
|
vars
|
-
|
(オプション) ODE 問題の未知数を示している、 1 変数の不定関数または 1 変数の不定関数の集合かリスト
|
procopts
|
-
|
(odesys が設定されていない場合に必要) プロシージャにより定義されたシステムを指定するためのオプション (procedure、initial、start、number、procvars)。詳細は、dsolve[numeric,IVP] を参照してください。
|
options
|
-
|
(オプション) keyword = value の形をした方程式
|
|
|
|
|
モデルの説明
|
|
|
初期および境界値の両問題に加え、微分代数の初期値 (DAE) 問題も数値的に解くことができます。システムが入力として与えられると、問題の種類は自動的に検出されますが、いくつかの DAE 問題では手法を指定する必要がある場合があります。
|
|
デフォルトでは、独立変数の値が与えられると、数値解を得るために使用できるプロシージャが返されます。
|
•
|
初期値問題 (IVP) に対するデフォルトの手法は、5 次精度の解を導き出す Runge-Kutta Fehlberg 法です(詳細については、rkf45 または numeric,IVP を参照)。境界値問題 (BVP) に対しては、 Richardson 外挿法を用いた有限差分法が使用されます (詳細については、numeric,BVP を参照)。微分代数の IVP 問題 (DAE) に対しては、rkf45 法の変形が使用されます (詳細については、dae_extension を参照)。
|
•
|
(デフォルトの stiff および nonstiff IVP 法については、 complex オプションを使用して複素問題であることを指定する必要がある場合がありますが) すべての IVP 法は、実数の独立変数を持つ複素数値の IVP を解くために使用できます。BVP または DAE 法については、現在、どちらも複素数値の問題には使用できません。そのため、dsolve を呼び出す前にシステムを実数値のシステムに変換しておく必要があります。
|
•
|
IVP および DAE 法は、返された解を求めるプロシージャに対して追加の機能を持っています。具体的には、最後に計算された解の値 (特異点を持つ問題に対して有用) の検索、初期データの検索および特定の解を求めるプロシージャについては初期データの変更 (DAE 法ではこの機能のサブセットをサポート) などができます。詳細については、numeric,IVP および numeric,DAE を参照してください。
|
•
|
従属変数の集合がリストとして指定された場合、リストにおける従属変数の順序が出力における従属変数の順序になります。それ以外の場合は、従属変数はその名前のアルファベット順に出力されます。
|
|
オプション
|
|
•
|
そのほとんどが IVP、 BVP および DAE に共通のハイレベルオプションは以下の通りです。
|
'output' =
|
keyword (キーワード) または array (配列)
|
'stiff' =
|
boolean (ブール値)
|
'events' =
|
list (リスト)
|
'event_pre' =
|
keyword (キーワード)
|
'event_maxiter' =
|
integer (整数値)
|
'range' =
|
numeric (数値)..numeric (数値)
|
'abserr' =
|
numeric (数値) または list (リスト)
|
'relerr' =
|
numeric (数値)
|
'known' =
|
name (名前) または list of names (名前のリスト)
|
'optimize' =
|
boolean (ブール値)
|
|
|
|
'output'= keyword または array
|
|
output は、procedurelist、listprocedure、operator または piecewise の値をとるキーワードか、数値解を求める点での独立変数の値を与える array または Array を値としてとります。
|
|
デフォルトのキーワードは、procedurelist で、dsolve からの出力をプロシージャとして返します。このプロシージャに引数として独立変数の値を与えると、解の値のリストが variable=value の形で戻ります。式で、左辺は独立変数、従属変数および (高階方程式の場合は) 導関数の名前を表し、右辺は対応する計算で得られた数値解を表します。
|
|
listprocedure キーワードは、出力を variable=procedure の形の等式のリストとして返します。式で、左辺は独立変数、従属変数および導関数を表し、右辺は対応する解の構成要素を計算する際に使用できるプロシージャを表します。この出力形態は、返されたプロシージャを後にfsolve 等で使用する場合などに最も有用です。
|
|
operator キーワードは、出力を operator=procedure の形の等式のリストとして返します。式で、左辺は、ある点で関数または導関数を評価するために与える独立変数の値に適用することができる演算子で、右辺は、 対応する解の構成要素の計算に使用することができるプロシージャを表します。この出力形態は、返されるリストをある点でショートカット評価 (shortcut evaluation) する際に最も有用です (たとえば、)。
|
|
piecewise キーワードは、non-stiff および stiff のデフォルトの IVP、DAE 法 (rkf45、ck45、rosenbrock、 rkf45_dae、ck45_dae、rosenbrock_dae) および taylorseries 法でのみ利用可能です。このキーワードは、出力を の形の等式のリストとして返します。式で、左辺は独立変数、従属変数および導関数を表し、右辺は対応する解の構成要素を記述する区分的多項式関数を表します。これらの区分的関数は、計算の各ステップに対する手法の補間式から得られています。さらに、区分的多項式の形は、区分的出力要求に対してインデックスを使用して指定することができます。output=piecewise[horner] を指定すると、出力は ホーナー (horner) 形式 (taylorseries (テイラー級数) に対するデフォルト) の多項式になります。一方、output=piecewise[polynomial] を指定すると、出力は標準の多項式形 (その他すべての手法に対するデフォルト) になります。この出力形態を使用する場合は、'range' 引数を用いて、希望する区分的関数の範囲を dsolve/numeric に指定する必要があります。
|
|
注意:Array が使用される場合は、出力に新しい Array と Matrix データタイプが用いられる以外、その出力は array の場合と同じです。出力点、あるいは解の成分が多く、インライン表示が Array または Matrix を表すプレースホルダーに変わるデフォルトのサイズを超えていると、出力データは直接表示されません。詳細については、interface(rtablesize) を参照してください。
|
|
'range'= numeric...numeric
|
|
- 解の変化が遅い広範囲な領域と、解の変化が急激な小さな領域を持つような問題について数値解を計算する場合には、'range' に odeplot の refine オプション (rkf45*、ck45*、rosenbrock* に限る) を組み合わせて使用することで、解の詳細をより良くプロットすることができます。'range' が使用されていない場合、dsolve の呼び出しは、計算した解を保存しないが、点が与えられたときに解を計算するプロシージャを返します。
|
|
- 長時間にわたる積分問題については、解全体の保存にかなりのメモリ容量が使用されるため 'range' を使用しないことを推奨します。
|
|
BVP 法に関しては、情報は deqns に含まれる境界条件から推測できるため、大域的な誤差限界の下で IVP を解く場合にのみこのオプションが必要になります。詳細については、numeric,BVP を参照してください。
|
|
'abserr'= numeric または list
|
|
IVP および DAE の場合には、成功したステップでの絶対誤差の許容限度を指定し、(例外的なシステムを除くすべての) BVP の場合には厳密解と計算値の最大誤差の目安を指定する数値です。この数値は、(ステップのサイズは固定で誤差制御がない) classical (古典) 法以外のすべての手法でサポートされています。IVP および DAE に対する rkf45*、ck45*、rosenbrock*、mebdfi 法では、各従属変数または各解成分に対して異なる絶対誤差の許容限度が指定できるリスト形式で abserr を利用することができます。詳細については、dsolve[numeric,IVP] および dsolve[Error_Control] のページを参照してください。
|
|
'abserr' および 'relerr' のデフォルトは、それぞれの手法に固有ですが、non-stiff の rkf45 および ck45 法については、デフォルトの許容誤差は (Digits の値で制御されるのではなく) 現在、固定であることに注意してください。
|
|
'known'= name または list of names
|
|
方程式システムを分析するときに既知と考える必要のある関数リストを提供します。これにより、ユーザ定義の関数を含む方程式システムの指定が可能になります。ユーザ定義の各関数は、数値入力が与えられた場合は常に数値を返し、名前の入力が与えられた場合は評価されずに戻る必要があります。このオプションの使用例については、下記アプリケーションと例題のセクションを参照してください。
|
|
注意:多くの場合、ユーザ定義の関数を含む方程式システムの使用は、dsolve 内でのevalhf の使用を妨げます。従って、この機能を使用して解を求めた場合、既知の関数に関してシステムを指定した場合 (この形でシステムを表現することが可能な場合) よりかなり多くの時間がかかることになります。
|
|
このオプションは、taylorseries、rosenbrock および rosenbrock_dae を除くすべての手法で直接サポートされていますが、これら 3 つの手法に関しては追加の情報を必要とします。主な問題は、これら 3 つの手法では入力 ODE システムの導関数の情報が必要で、関数がプロシージャとして定義されている場合はこれらを直接利用できないことです。ただし、関数の導関数を 1 つまたは複数の `diff/` ルールを用いて定義し、'known' を指定することによりこれらの手法が使用可能になります。rosenbrock 法は 1 つの導関数のみ必要としますが、taylorseries 法は任意階の導関数を必要とします。そのため、taylorseries はプロシージャではなく、関数として表現可能で、その関数の導関数に至るまですべての導関数に対して diff ルールが定義されている 'known' 関数の有限階導関数が存在する場合にのみ有用です。rosenbrock_dae 法は、DAE の前処理で入力システムを微分する必要がある場合があるため、1 つ以上の導関数を必要とする場合があります。
|
|
dsolve に、入力システムを計算のために最適化するよう命令するブール値です。最適化を行うと、処理は増えますが、計算の高速化が期待できます。 このオプションのデフォルト値は、IVP および BVP 問題に対しては、false で、(通常、一番効果が上がる) DAE 問題に対しては、true です。 計算に要する時間を長引かせる場合があるため、(piecewise などの) 条件付き評価関数を含むシステムでは optimize=true は使用しないでください。このオプションは、システム定義の問題および (独自の最適化を使用する) taylorseries を除くすべての手法で利用可能です。
|
|
|
注意
|
|
•
|
Digits<= の場合、数値解は倍精度 (機械精度の浮動小数) で計算されます。Maple の初期設定で、Digits のデフォルトは 10 に設定されています。明示的に、Digits> と設定した場合、dsolve は、機械精度の浮動小数の代わりに、Maple の浮動小数で数値解を計算するようになります。dsolve の呼び出しにおいて、計算の精度は固定です。そのため、出力として返されるプロシージャは、一旦作成されると、Digits の設定が変更されても計算された解に影響を与えることがありません。
|
|
数値解の計算には、fsolve を併用することができますが、その場合、fsolve の呼び出しにおける Digits は dsolve プロシージャの精度より (通常は、5 またはそれ以下) 低い精度に設定される必要があります。 (これは、fsolve は残りに対して現行の Digits 設定よりも高い精度を必要とするためです) rkf45、ck45、rosenbrock および taylorseries 法では、補間式で数値解を計算する場合は、fsolve から呼び出されたときの現行の Digits 設定を用いて補間式を計算することで問題を解決します。この問題の例については、下記最後の例を参照してください。
|
|
|
|
アプリケーションと例題
|
|
リストとして指定されている従属変数の集合:
>
|
dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0}: #default alphabetical
|
>
|
dsn1 := dsolve(dsys,numeric):
|
| (4.1) |
>
|
[t = 1., x(t) = 0.540302304342753614, y(t) = -0.841471136025939148];
|
| (4.2) |
>
|
dsn2 := dsolve(dsys,numeric,[y(t),x(t)]): # specify differently
|
| (4.3) |
>
|
[t = 1., y(t) = -0.841471136025939148, x(t) = 0.540302304342753614];
|
| (4.4) |
単純な 2 階 IVP:
>
|
deq1 := (t+1)^2*diff(y(t),t,t) + (t+1)*diff(y(t),t)
+ ((t+1)^2-0.25)*y(t) = 0;
|
| (4.5) |
>
|
ic1 := y(0) = 0.6713967071418030,
D(y)(0) = 0.09540051444747446:
|
>
|
dsol1 := dsolve({deq1,ic1}, numeric, range=0..1);
|
| (4.6) |
| (4.7) |
| (4.8) |
| (4.9) |
出力が演算子の同じ 2 階 IVP:
>
|
dsol1b := dsolve({deq1,ic1}, numeric, output=operator);
|
| (4.10) |
| (4.11) |
| (4.12) |
| (4.13) |
1 階 IVP:
>
|
dsol2 := dsolve({diff(y(x),x) = y(x)*cos(x), y(0) = 1},
numeric, output=listprocedure);
|
| (4.14) |
>
|
fy2 := eval(y(x),dsol2);
|
| (4.15) |
>
|
[fy2(0), fy2(Pi/4), fy2(Pi/2), fy2(3*Pi/2), fy2(Pi)];
|
| (4.16) |
区分出力を伴う 1 階 IVP:
>
|
dsol2p := dsolve({diff(y(x),x) = y(x)*cos(x), y(0) = 1},
numeric, output=piecewise, range=0..1);
|
| (4.17) |
硬い (stiff) 非線形問題:
>
|
dsys3 := diff(x(t),t) = x(t)^2 - t^6, x(0)=-0.1;
|
| (4.18) |
>
|
dsol3 := dsolve([dsys3], type=numeric, stiff=true,
output=Array([0,0.25,0.5,0.75,1]));
|
| (4.19) |
線形 BVP:
>
|
dsol4 := dsolve([diff(y(x),x,x) = 3*y(x),
y(0) = 1, y(2) = 1.2], numeric);
|
| (4.20) |
| (4.21) |
| (4.22) |
| (4.23) |
| (4.24) |
非線形 BVP:
>
|
dsys5 := {diff(y(x),x,x) + abs(y(x)) = 0, D(y)(0) = 1, y(1) = -1};
|
| (4.25) |
>
|
dsol5 := dsolve(dsys5, numeric, output=Array([0,0.25,0.5,0.75,1]));
|
| (4.26) |
DAE:
>
|
dsys6:= {
diff(x(t),t,t) = -2*lambda(t)*x(t),
diff(y(t),t,t) = -2*lambda(t)*y(t)-Pi^2,
x(t)^2+y(t)^2 = 1,
x(0)=0, D(x)(0)=1/10, y(0)=-1, D(y)(0)=0};
|
| (4.27) |
>
|
dsol6 := dsolve(dsys6, numeric);
|
| (4.28) |
| (4.29) |
fsolve の使用:
>
|
dsys7 := {diff(y(x),x)=2*y(x), y(0)=1};
|
| (4.30) |
| (4.31) |
>
|
dsol7 := dsolve(dsys7, numeric, method=gear, abserr=1e-16, relerr=1e-16,
output=listprocedure);
|
| (4.32) |
>
|
yproc := rhs(dsol7[2]);
|
| (4.33) |
同じ Digits 設定での試行 - 評価なしでの戻り
>
|
fsolve(yproc(x)=1e5, x=4..7);
|
| (4.34) |
Digits を 15 に設定した上での同じコマンドの実行 - 成功
| (4.35) |
>
|
fsolve(yproc(x)=1e5, x=4..7);
|
| (4.36) |
| (4.37) |
known オプションの使用:
>
|
f := proc(x) local t;
if not type(evalf(x),'numeric') then
'procname'(x);
else
evalf(Int(exp(-t^2/10),t=0..x));
end if;
end proc;
|
| (4.38) |
入力が数値でない場合、プロシージャは評価せずに戻るよう設定されていることに注意してください。
>
|
dsys8 := {diff(y(x),x)=y(x)+f(x), y(0)=0};
|
| (4.39) |
>
|
dsolve(dsys8, numeric);
|
エラー、(`dsolve/numeric/DAE/make_proc` で) 未知関数と方程式の数は一致している必要がありますが、関数が 2 つ {f, y} と方程式が 1 つしか指定されていません。
| |
'known' が指定されていなかったため、失敗しました。次では、指定されています。
>
|
dsol8 := dsolve(dsys8, numeric, known=f);
|
| (4.40) |
| (4.41) |
| (4.42) |
導関数を必要とする手法:
>
|
dsolve(dsys8, numeric, method=rosenbrock, known=f);
|
エラー、(`dsolve/numeric/checkknown` で) プロシージャ定義の関数 f で 'rosenbrock' 法を使用するには、微分のルール `diff/f` が必要です。
| |
微分のルールが必要であるため、失敗しました。diff の記憶テーブルを消去し、次のように定義します:
>
|
`diff/f` := proc()
# Chain rule
diff(args[1],args[2])*exp(-(args[1])^2/10);
end proc:
|
これは、次のように動作します:
| (4.43) |
>
|
diff(Int(exp(-t^2/10),t=0..sin(x)),x);
|
| (4.44) |
今度は、次のようになります:
>
|
dsol8b := dsolve(dsys8, numeric, method=rosenbrock, known=f);
|
| (4.45) |
| (4.46) |
| (4.47) |
|
|
関連項目
|
|
Digits, dsolve/dae_extension, dsolve/Error_Control, dsolve[ck45], dsolve[classical], dsolve[dverk78], dsolve[Events], dsolve[gear], dsolve[interactive], dsolve[lsode], dsolve[maxfun], dsolve[mebdfi], dsolve[numeric,BVP], dsolve[numeric,DAE], dsolve[numeric,interactive], dsolve[numeric,IVP], dsolve[rkf45], dsolve[rosenbrock], dsolve[Stiffness], dsolve[taylorseries], fsolve, interface, plots[odeplot]
|
|
リファレンス
|
|
|
Ascher, U.; Mattheij, R.; and Russell, R. "Numerical Solution of Boundary Value Problems for Ordinary Differential Equations." SIAM Classics in Applied Mathematics. Vol. 13. (1995).
|
|
Ascher, U., and Petzold, L. "Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations." SIAM, Philadelphia. 1998.
|
|
Bailey, P.B.; Shampine, L.F.; and Waltman, P.E. Nonlinear Two Point Boundary Value Problems. New York: Academic Press, 1968.
|
|
Boyce, W.E., and DiPrima, R.C. Elementary Differential Equations and Boundary Value Problems. New York: John Wiley & Sons, 1997.
|
|
Cash, J.R. "The Integration of Stiff IVP in ODE Using Modified Extended BDF." Computers and Mathematics with Applications. Vol. 9. (1983): 645-657.
|
|
Gear, C.W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, 1971.
|
|
Hindmarsh, Alan C.; Stepleman, R.S.; et al, eds. Odepack, a Systemized Collection of ODE Solvers. Amsterdam: North-Holland, 1983.
|
|
Hubbard J.H., and West, B.H. Differential Equations: A Dynamical Systems Approach. New York: Springer, 1990. Part I. One Dimensional Equations.
|
|
Hull, T.E.; Enright, W.H.; Fellen, B.M.; and Sedgwick, A.E. "Comparing Numerical Methods for Ordinary Differential Equations." SIAM J. Numer. Anal. Vol. 9. (1972): 603-637.
|
|
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.
|
|
|