|
dsolve/numeric/ck45 - 常微分方程式 (ODE)の数値解を求める
|
|
使い方
|
|
dsolve(odesys, numeric, method=ck45, vars, options)
dsolve(numeric, method=ck45, procopts, options)
|
|
|
パラメータ
|
|
|
odesys
|
-
|
集合またはリスト；常微分方程式および初期条件
|
|
numeric
|
-
|
名前；dsolve に数値解を要求します
|
|
method=ck45
|
-
|
方程式：数値解法として使用します
|
|
vars
|
-
|
(オプション)常微分方程式の未知変数を含む全ての一変数の不定関数, あるいはそれらの集合かリスト
|
|
options
|
-
|
(オプション) keyword = value 形式の方程式
|
|
procopts
|
-
|
常微分方程式に用いる手続きの指定オプション (procedure, initial, start, number, procvars のいずれか)。詳細は dsolve/numeric/IVP をご参照ください。
|
|
|
|
|
|
モデルの説明
|
|
|
•
|
オプションに numeric および method=ck45 を指定した dsolve コマンドは、4次の補間を含む 4-5 次 Runge-Kutta Cash-Karp 法(Cash-Karp fourth-fifth order Runge-Kutta)を用いて数値解を求めます。stiff を使用しない場合、この解法は初期値問題(IVP)を type=numeric と指定して解く場合の、デフォルトが dsolve/rkf45 であるもうひとつの non-stiff の解法です。
|
|
|
操作方法
|
|
ck45 による解法には、二つの異なった(手続き型出力のための)操作方法があります。
|
|
range オプションを用いる場合、この解法は、指定範囲(range)において、返された手続きを全て呼び出すために要求された解の値を高速に補間し、それを解の情報として内部に記憶して初期値問題(IVP)の解を計算します。
|
|
|
指定範囲外の点で呼び出された手続きを返すことも可能ですが、あまり薦められません。
|
|
|
この方法は、適応型プロットとして odeplot の refine オプションと組み合わせて使用することができます(つまり、refine が指定されたとき、odeplot はプロットのために前計算された点を使います) 。
|
|
|
この方法は、解が特異になれるような問題に使用することは薦められません。各ステップを保存すると、特異性の近くでは多くのステップを取ることになり、メモリの使用量が深刻な問題となります。
|
|
|
この方法を使用している間に記憶する補間の情報は、以下で説明するようように interpolation=false オプションを用いることで無効にすることができます。これは、(離散解に加えて)補間の記憶により多くのメモリを必要とする場所で高精度な解を得るために薦められます。解の値は最も近い 5 点の補間から得られるため、補間を無効にすることは一般的に薦められません。さらに、必ずしも 4 つのエラーを補間に提供するわけではありません。
|
|
|
range オプションを用いない場合、初期値問題(IVP)の解の値は記憶されませんが、要求される度に計算することになります。
|
|
|
全ての解の値が記憶されるわけではないので、初期点と最も新しく計算された点の間の点が要求された場合はいつも、( 積分の反転を避けるために ) 初期値から計算しなおさなければなりません。したがって、初期値から移動した解の値を集めることが賢明です。
|
|
|
|
オプション
|
|
|
•
|
T以下のオプションは ck45 で使用することが可能です。
|
|
'output'
|
=
|
keyword or array
|
|
'known'
|
=
|
name or list of names
|
|
'abserr'
|
=
|
numeric
|
|
'relerr'
|
=
|
numeric
|
|
'initstep'
|
=
|
numeric
|
|
'interr'
|
=
|
boolean
|
|
'maxfun'
|
=
|
integer
|
|
'number'
|
=
|
integer
|
|
'procedure'
|
=
|
procedure
|
|
'start'
|
=
|
numeric
|
|
'initial'
|
=
|
array
|
|
'procvars'
|
=
|
list
|
|
'startinit'
|
=
|
boolean
|
|
'implicit'
|
=
|
boolean
|
|
'optimize'
|
=
|
boolean
|
|
'compile'
|
=
|
boolean or auto
|
|
'range'
|
=
|
numeric..numeric
|
|
'events'
|
=
|
list
|
|
'event_pre'
|
=
|
keyword
|
|
'event_maxiter'
|
=
|
integer
|
|
'event_iterate'
|
=
|
keyword
|
|
'event_initial'
|
=
|
boolean
|
|
'complex'
|
=
|
boolean
|
|
|
|
|
dsolve の出力形式を指定します。procedurelist, listprocedure, または operator は手続き型が出力を指定し、piecewise は指定された範囲の独立変数の値で piecewise 関数の形式で出力することを指定します。さらに 1-D 配列および配列は固定された独立変数の値を出力します。詳細は dsolve/numeric をご参照ください。
|
|
|
abserr, relerr および initstep
|
|
|
解の精度と、解法の最初のステップ幅を指定します。詳細は Error_Control をご参照ください。ck45 でのデフォルトの値は abserr=1e-7 と relerr=1e-6 となります。もし指定しなかった場合、initstep の値は連立常微分方程式の局所的な振る舞いを考慮して、解法が決定します。
|
|
|
この値は、デフォルトでは true に設定されており、(常微分方程式の index-1 変数の補間を含む)解の補間誤差を誤差管理に統合するかどうかを制御します。 に設定すると、(piecewise による導関数の不連続点など)解が急速に変化する場所には、指定された許容誤差により指示されるものより遥かに大きな誤差が発生する可能性があります。 に設定すると、これらの範囲における誤差を最小化するためにステップ幅は縮小されますが、変数にジャンプ不連続点がある問題では、積分は特異性があることを示す誤差の範囲に入ります。誤差が発生する後者の場合、events (dsolve/Events 参照) を使用した不連続点のモデル化の方が利点があるかもしれません。
|
|
|
一階の連立常微分方程式の右辺の評価数の最大値を指定します。このオプションは  と指定することで無効にできます。ck45 でのデフォルトの値は  です。
|
|
|
number, procedure, start, initial および procvars
|
|
|
startinit, implicit, および optimize
|
|
|
このオプションは、数値解の計算に使用される内部で作成された手続きをコンパイルして効率化するよう指定します。このオプションは、 がハードウェアの精度範囲内に設定され、入力関数が evalhf が可能な関数のみ(例えば、exp, sin, および ln のような初歩的数学関数のみ)を含んでいる場合に限って、機能することに注意してください。dsolve/Efficiency をご参照ください。デフォルトではこの値は false に設定されています。この値が true に設定されていてコンパイルが不可能な場合、エラーが表示されたます。auto に設定されていてコンパイルが不可能な場合は、コンパイルされない手続きが直接使用されます。
|
|
|
解の値が必要とする独立変数の値の範囲を決定します。このオプションの使用は、dsolve/numeric に説明がある手続き型の出力形式 (output types) のために計算法の振る舞いを大きく変えることになります (このページの操作方法 ( Modes of Operation ) の range を使用する場合としない場合についての説明をご参照ください ) 。
|
|
|
events, event_pre, event_maxiter, event_iterate, event_initial
|
|
|
これらのオプションを使用して、数値解のイベントの振る舞いを指定および制御します。これらのオプションについては、dsolve/Events で詳細に説明されています。
|
|
|
問題が複素数値である(またはそうなるような)場合を表示する論理演算子を受け入れます。デフォルトでは方程式系と初期データの入力に基づいて検出するように設定されていますが、入力した方程式系が手続き型で定義されている場合や、初期データが実数である場合には、解を得るために complex=true を指定する必要があります。これは最初の実方程式系が複素方程式系になる可能性があります。この変化が起こる点は特異であると考えられるので、もし complex=true が指定されていなければ、この点で積分が停止することになります。
|
|
|
|
追加項目
|
|
|
•
|
infolevel[dsolve] を  と設定することで、最後の格子点での評価値の情報に誤差が含まれます。
|
|
•
|
ck45 による方法は、Digits の環境変数を変更して計算精度を上げることで、連立初期値問題 (IVPs) の解を高精度に計算することが可能になります。しかしながら、これは高階な dverk78 か gear を使う場合に薦められる方法です。
|
|
•
|
結果は plots パッケージの odeplot 関数を使用してプロットすることができます。
|
|
|
|
|
アプリケーションと例題
|
|
範囲を指定した解 - 区間  の解を記憶します。：
|
>
|
dsys1 := {diff(x(t),t)=y(t), diff(y(t),t)=x(t)+y(t),
x(0)=2, y(0)=1};
|
|
| (4.1) |
|
>
|
dsol1 := dsolve(dsys1, numeric, method=ck45, output=listprocedure,
range=0..1):
|
|
>
|
dsol1x := subs(dsol1,x(t)): dsol1y := subs(dsol1,y(t)):
|
|
>
|
[dsol1x(0), dsol1y(0)];
|
|
![[2., 1.]](/support/helpjp/helpview.aspx?si=1995/file00995/math863.png)
| (4.2) |
|
>
|
[dsol1x(0.4), dsol1y(0.4)];
|
|
![[2.69118406139209, 2.60811663032515]](/support/helpjp/helpview.aspx?si=1995/file00995/math870.png)
| (4.3) |
|
>
|
[dsol1x(1.0), dsol1y(1.0)];
|
|
![[5.58216687757548, 7.82688820381604]](/support/helpjp/helpview.aspx?si=1995/file00995/math877.png)
| (4.4) |
範囲を指定しない解 - 各点では要求されたときに計算します。：
|
>
|
dsol2 := dsolve(dsys1, numeric, method=ck45,
output=procedurelist):
|
|
![[t = 0., x(t) = 2., y(t) = 1.]](/support/helpjp/helpview.aspx?si=1995/file00995/math893.png)
| (4.5) |
|
![[t = .4, x(t) = 2.69118406139209, y(t) = 2.60811663032515]](/support/helpjp/helpview.aspx?si=1995/file00995/math900.png)
| (4.6) |
|
![[t = 1.0, x(t) = 5.58216687757548, y(t) = 7.82688820381604]](/support/helpjp/helpview.aspx?si=1995/file00995/math907.png)
| (4.7) |
配列の解：
|
>
|
dsys3 := {(D@@2)(x)(t)=-y(t), (D@@2)(y)(t)=D(x)(t)+y(t)};
|
|
| (4.8) |
|
>
|
init3 := {x(0)=1, D(x)(0)=0, y(0)=0, D(y)(0)=1};
|
|
| (4.9) |
|
>
|
dsol3 := dsolve(dsys3 union init3, numeric, method=ck45,
output=Array([0,.6,1.1,1.5,2.3,2.5]));
|
|
| (4.10) |
procedurelist 出力：
|
>
|
deqn4 := { diff(y(t),t$3) - 2*diff(y(t),t$2) + 2*y(t) };
|
|
| (4.11) |
|
>
|
init4 := { y(0)=1, D(y)(0)=1, (D@@2)(y)(0)=1 };
|
|
| (4.12) |
|
>
|
dsol4 := dsolve(deqn4 union init4, numeric, method=ck45,
relerr=Float(1,-8), abserr=Float(1,-8),
maxfun=0, output=procedurelist);
|
|
| (4.13) |
|
![[t = 0.5e-1, y(t) = 1.05124946387462, diff(y(t), t) = 1.04995674032366, diff(y(t), `$`(t, 2)) = .997371824249566]](/support/helpjp/helpview.aspx?si=1995/file00995/math966.png)
| (4.14) |
|
![[t = -.45, y(t) = .648630780960807, diff(y(t), t) = .571757531626938, diff(y(t), `$`(t, 2)) = .870897129164657]](/support/helpjp/helpview.aspx?si=1995/file00995/math973.png)
| (4.15) |
複素数値問題における Piecewise 出力：
|
>
|
deqn5 := { diff(y(t),t) = I*y(t) };
|
|
| (4.16) |
|
| (4.17) |
|
>
|
dsol5 := dsolve(deqn5 union init5, numeric, method=ck45,
output=piecewise, complex=true, range=0..Pi):
|
|
| (4.18) |
|
| (4.19) |
|
>
|
eval(dsol5[2],t=evalf[15](Pi));
|
|
| (4.20) |
|
|
|
関連項目
|
|
dsolve/classical, dsolve/dverk78, dsolve/Error_Control, dsolve/Events, dsolve/gear, dsolve/lsode, dsolve/maxfun, dsolve/numeric, dsolve/numeric/IVP, dsolve/rkf45, dsolve/rosenbrock, dsolve/Stiffness, dsolve/taylorseries, infolevel, plots[odeplot]
|
|
|
リファレンス
|
|
|
|
Enright, W.H.; Jackson, K.R.; Norsett, S.P.; and Thomsen, P.G. "Interpolants for Runge-Kutta Formulas." ACM TOMS, Vol. 12, (1986): 193-218.
|
|
|
Cash, J.R.; Karp, A.H. "A Variable Order Runge-Kutta Method for 初期値 Problems with Rapidly Varying Right-Hand Sides". ACM Transactions on Mathematical Software, Vol. 16(3), (1990): 201-222.
|
|
|
