dsolve/numeric/classical - 常微分方程式の数値解を求める
使い方
dsolve(odesys, numeric, method=classical)
dsolve(odesys, numeric, method=classical[choice], vars, options)
dsolve(numeric, method=classical[choice], procopts, options)
パラメータ
odesys - 集合またはリスト ; ( 連立 ) 常微分方程式および初期条件
numeric - 名前 ; dsolve に数値解を要求します
method=classical - method=classical を指定します ; 数値解法として使用します
method=classical[choice] - 方程式 ; 数値解法と共に使用する副解法の指定
vars - ( オプション ) odesys に現れる従変数、または従変数の集合かリスト
options - ( オプション ) keyword = value 形式の方程式
procopts - 連立常微分方程式に用いる手続きの指定オプション (procedure, initial, start, number, procvars のいずれか。 ) 。詳細は dsolve[numeric,IVP] をご参照ください。
|
説明
|
|
•
|
オプションに numeric および method=classical または method=classical[choice] を指定した dsolve コマンドは、以下で説明されるような古典的な数値解法 (classical numerical methods) を用いて数値解を求めます。
|
|
この計算法はステップ幅を固定し、許容誤差および修正誤差を全く与えずに使用します。これは主に教育で扱うために提供します。もし、常微分方程式の初期値問題 (IVP) の数値解法の練習として必要ならば、dsolve[numeric,IVP] のヘルプページで紹介した全ての計算法の代わりに使用するように指示してください。
|
•
|
method=classical[choice] を選択した古典的計算法は choice オプションの使用を通して利用可能です。デフォルトでは foreuler となっており、これは前進 Euler 法 (forward Euler method) のことを意味します。
|
•
|
利用可能な古典的 (classical) 計算法の選択は、y'=f(t,y) なる問題に対して説明されます。ここで Y[i] は時刻 t[i] での解の許容値であり、h は固定されたステップ幅 t[i]-t[i-1] で、さらに時刻 t[n+1] での解 Y[n+1] の値は各公式によって計算されます。
|
•
|
foreuler は方程式で定義された前進 Euler 法 :
|
Y[n+1] = Y[n] + h*f(t[n], Y[n])
•
|
heunform は Heun 公式 (Heun formula : 改良 Euler 法 (improved Euler method) として知られています。 ) を意味します。これは解を予測するために前進 Euler 法を使用し、修正として台形則 (trapezoid rule) を受け入れます。これらは方程式のペアとして指定されています。:
|
Yp = Y[n] + h*f(t[n], Y[n])
Y[n+1] = Y[n] + (h/2)*(f(t[n],Y[n]) + f(t[n+1],Yp))
•
|
impoly は修正多角形法 (improved polygon method : 修正 Euler 法 (modified Euler method) として知られています。) を意味します。これは方程式として指定されています。:
|
Y[n+1] = Y[n] + h*(f(t[n]+h/2, Y[n]+(h/2)*f(t[n],Y[n])))
•
|
rk2 は二次の古典的 Runge-Kutta 法 (second-order classical Runge-Kutta method) を意味します。これは以下で指定されています。:
|
k1 = f(t[n], Y[n])
k2 = f(t[n]+h, Y[n]+h*k1)
Y[n+1] = Y[n] + (h/2)*(k1+k2)
|
注意 : これは Heun 公式と同じものとなります。
|
•
|
rk3 は三次の古典的 Runge-Kutta 法 (third-order classical Runge-Kutta method) を意味します。これは以下で指定されています。:
|
k1 = f(t[n], Y[n])
k2 = f(t[n]+(h/2), Y[n]+(h/2)*k1)
k3 = f(t[n]+h, Y[n]+h*(-k1+2*k2))
Y[n+1] = Y[n] + (h/6)*(k1+4*k2+k3)
•
|
rk4 は四次の古典的 Runge-Kutta 法 (fourth-order classical Runge-Kutta) を意味します。これは以下で指定されています。:
|
k1 = f(t[n], Y[n])
k2 = f(t[n]+h/2, Y[n]+(h/2)*k1)
k3 = f(t[n]+h/2, Y[n]+(h/2)*k2)
k4 = f(t[n]+h, Y[n]+h*k3)
Y[n+1] = Y[n] + (h/6)*(k1+2*k2+2*k3+k4)
|
これは method=rkf45 による計算法と混乱しないように注意してください。method=rkf45 は 4-5 次 Runge-Kutta Fehlberg 法 (Fehlberg fourth-fifth order Runge-Kutta method) を使用しています。
|
•
|
adambash は Adams-Bashforth 法 (Adams-Bashforth method : 予測子法 (predictor method) として知られています。) を意味しています。これは以下で指定されています。:
|
Y[n+1] = Y[n] + (h/24) * (55*f(t[n],Y[n]) - 59*f(t[n-1],Y[n-1])
+ 37*f(t[n-2],Y[n-2]) - 9*f(t[n-3],Y[n-3]))
•
|
abmoulton は Adams-Bashforth-Moulton 法 (Adams-Bashforth-Moulton method : 予測子・修正子法 (predictor-corrector method) として知られています。) を意味します。これは以下で指定されています。:
|
Y[n+1] = Y[n] + (h/24) * (9*f(t[n+1],Y[n+1]) + 19*f(t[n],Y[n])
- 5*f(t[n-1],Y[n-1]) + f(t[n-2],Y[n-2])),
|
最初に Adams-Bashforth 法 ( 予測子法 ) によって f(t[n+1],Y[n+1]) が発見され、その後に Adams-Bashforth-Moulton 法 ( 予測子法 ) を使用します。
|
•
|
以下のオプションは classical で使用することが可能です。:
|
'output' = keyword or array
'known' = name or list of names
'maxfun' = integer
'number' = integer
'procedure' = procedure
'start' = numeric
'initial' = array
'procvars' = list
'startinit' = boolean
'implicit' = boolean
'optimize' = boolean
'stepsize' = numeric
'corrections' = integer
|
dsolve の出力形式を指定します。さらに、known オプションはユーザが定義した既知の関数を指定します。詳細は dsolve[numeric] をご参照ください。
|
|
一階の連立常微分方程式の右辺の評価数の最大値を指定します。このオプションは maxfun=0 と指定することで無効にできます。classical でのデフォルトの値は 50000 です。詳細は dsolve[maxfun] をご参照ください。
|
|
'number', 'procedure', 'start', 'initial', and 'procvars'
|
|
'startinit','implicit', and 'optimize'
|
|
各ステップの静的な幅を指定します。注意として、classical ではステップ幅の適合による許容誤差の修正を使用しません。デフォルトのステップ幅は、前進積分に対しては h = min((T0-Tf)/3, 0.005) となっており、後退積分に対しては h = max((T0-Tf)/3, -0.005) となっています。それぞれの問題は、T0 から Tf.までの間隔で解かれています。
|
|
各ステップで abmoulton 法によって受け入れられた修正点の数を正の整数として指定します。これは 4を超えないように設定することを薦めます ( 生成子の値が 4 以上である場合はより性格な結果を与えることができず、さらに長い時間を費やしてしまうことになります。) 。
|
|
|
例
|
|
>
|
dsys1 := diff(y(t),t$3) = y(t)+diff(x(t),t),
diff(x(t),t$2) = x(t)*y(t)-diff(y(t), t$2);
|
| (2.1) |
>
|
init1 := y(0) = 1, D(y)(0) = 2, (D@@2)(y)(0) = -1,
x(0) = 4, D(x)(0) = 4.3;
|
| (2.2) |
>
|
ans1 := dsolve({dsys1,init1}, numeric,
method=classical[abmoulton], corrections=2);
|
| (2.3) |
| (2.4) |
>
|
deq2 := diff(y(x), x$2) = y(x);
|
| (2.5) |
>
|
init2 := y(0) = 1, D(y)(0) = 1;
|
| (2.6) |
>
|
Digits := 20:
ans2 := dsolve({deq2, init2}, numeric,
method=classical[heunform],
output=array([0.01,0.1,1]),
stepsize=0.001);
|
| (2.7) |
>
|
Digits := 10:
deq3 := { diff(y(t),t$4) - diff(y(t),t$3) = y(t)*t^2 };
|
| (2.8) |
>
|
init3 := { y(0) = 3.56, D(y)(0) = 12, (D@@2)(y)(0) = -4,
(D@@3)(y)(0) = 6.544 };
|
| (2.9) |
>
|
ans3 := dsolve(deq3 union init3, numeric,
method=classical[rk3], output=listprocedure):
sol3 := eval(y(t), ans3);
|
| (2.10) |
| (2.11) |
| (2.12) |
|
|
参照
|
|
dsolve[dverk78], dsolve[gear], dsolve[lsode], dsolve[maxfun], dsolve[numeric], dsolve[numeric,IVP], plots[odeplot], dsolve[rkf45], dsolve[rosenbrock], dsolve[taylorseries]
|
|
参考文献
|
|
|
Boyce, W.E., and DiPrima, R.C. Elementary Differential Equations and Boundary Value Problems. 5th ed., 1992.
|
|
Conte, S.D., and C. de Boor. Elementary Numerical Analysis, An Algorithmic Approach. 1980.
|
|
Fox, L., and Mayers, D.F. Numerical Solution of Ordinary Differential Equations for Scientists and Engineers. 1987.
|
|
Lambert, J.D. Computational Methods in Ordinary Differential Equations. 1973.
|
|
|