derkf derkf_int

derkf() def derkf ( info  , n  , f  , t  , y  , tout  , wsave  , iwsave  , rtol  = 1.0e-10, atol  = 1.0e-10, mode  = 0  ) 常微分方程式の初期値問題 (5(4)次 ルンゲ・クッタ・フェールベルグ法) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは, 微分方程式ソルバー・パッケージDEPACに収録されているプログラムのひとつである(DEPACは, derkf, deabmおよびdebdfからなる). derkfは5(4)次ルンゲ・クッタ・フェールベルグ法のプログラムである. 3つの中ではアルゴリズムおよび使用法の両面で最もシンプルなルーチンである. 主に, 微分係数の計算に時間がかからない非スティフおよびややスティフな微分方程式を解くために設計されている. 一般的には, 高精度な解を求める場合や非常に多くの点における解を求める場合には適さない. derkfはオーバーヘッドが非常に小さいため, 中程度の精度が必要で方程式の計算に時間がかからない場合に, 通常は最も短い計算時間で解を求める. derkfは H. A. Watts および L. F. Shampine によるプログラム RKF45 の修正版のドライバー・ルーチンである. 下記参考文献を基にした密出力機能が提供される(modeおよび使用例(2)を参照せよ). 戻り値 (t1, info1) t1 (float): 独立変数tの最終ステップの最後の点の値 (インターバル・モードでは通常toutに等しい). 解がこの点まで正常に求められたことを示す. info1 (int): リターンコード. この値をチェックし, info1 = 1〜7 であれば必要により次のアクションとして t = t1, info = info1 と設定したうえで本ルーチンを再度呼び出すことができる. = -1: 入力パラメータ info の誤り (info < 0 または info > 7) = -2: 入力パラメータ n の誤り (n < 1) = -3: 入力パラメータ f の誤り = -4: 入力パラメータ t の誤り (継続呼び出しで t = tout または t != 前回のtout) = -5: 入力パラメータ y の誤り = -6: 入力パラメータ tout の誤り (積分の方向を変更した) = -7: 入力パラメータ wsave の誤り = -8: 入力パラメータ iwsave の誤り = -9: 入力パラメータ rtol の誤り (rtol < 0) = -10: 入力パラメータ atol の誤り (atol < 0) = 1: 正常終了 (t = tout). toutを再設定して再呼び出しすることができる. toutはtと異なる値であること. = 2: 中間結果出力モードのため戻った (まだ toutに達していない). tout方向の次のステップを続行するために再呼び出しすることができる. = 3: 最大ステップ数(10000)を超えた. 再呼び出しして続行することができる. さらに10000ステップ続けることができる. = 4: 誤差許容値が厳しすぎる. 許容値を変更してから再呼び出しすることができる. = 5: 求められた解が0であるため, 純粋な相対誤差テスト(atol = 0)ができない. atolを正の値にして再呼び出しして続行することができる. = 6: スティフな問題の可能性がある. 再呼び出しして続行することができる. しかし, スティフな問題を効率的に扱うことができるDEPACのdebdfのような他のプログラムが推奨される. = 7: 頻繁な出力のためステップサイズが制限された. 再呼び出しして続行することができる. しかし, 密出力機能を持った他のプログラムが推奨される. = 8: 無限ループが検出された. 引数 [in] info 制御コード. = 0: 最初の呼び出し時(新たに計算を開始する場合)には info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される. = 1〜7: info1 = 1〜7 で戻った場合, 計算を継続するために, info = info1 と設定して再呼び出しすることができる. [in] n 微分方程式の数. (n >= 1) [in] f 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. _CODE def f(n, t, y, yp): yp[i] = tおよびyにおける微分係数の計算値 (i = 0 〜 n-1) _ENDCODE ただし, nは方程式の数, yp(長さnの1次元配列)は与えられたtおよびy(長さnの1次元配列)における微分係数, すなわち, yp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 〜 n-1). [in] t 独立変数tの初期値. 本ルーチンはtからtoutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点がt1に返される. 新たなtoutにおける解を求めるために t = t1 として積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, tは前回のtoutに等しくなければならない. また, toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい. モードはパラメータmodeで指定できる. [in,out] y Numpy ndarray (1次元配列, float, 長さn) [in] tの初期値における従属変数yの初期値. [out] 最終のt(インターバル・モードにおいてはtoutに等しい)において求められた解(の近似値). [in] tout 解を求めたい点を設定する. tout = t とすることができ, その場合はtにおける微分係数を計算して戻る. 積分を行うのはtについて前進方向(tout > t)でも後退方向(tout < t)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてtからtoutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても従来どおりtoutを指定しなければならない. [out] wsave Numpy ndarray (1次元配列, float, 長さ 7*n + 20 (ただし, mode = 2 のときは 9*n + 20)) データ領域 (info = 0 として呼び出してから一連の計算の間変更してはならない). [out] iwsave Numpy ndarray (1次元配列, int, 長さ 20) 整数データ領域 (info = 0 として呼び出してから一連の計算の間変更してはならない). [in] rtol (省略可) 求める解の精度を指定する相対誤差許容値. (省略時 = 1.0e-10) 許容値は各ステップにおける局所誤差テストに使われ, yの各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).   abs(局所誤差) <= rtol*abs(y[i]) + atol rtol = 0 と設定すると純粋に絶対誤差テストとなる. 誤差が約1.0e-8以下の相対精度を必要とするならば, 通常はderkfを使わないほうがよい. DEPACのdeabmルーチンは厳しい精度を必要とする場合により効率的である. [in] atol (省略可) 求める解の精度を指定する絶対誤差許容値. (省略時 = 1.0e-10) 許容値は各ステップにおける局所誤差テストに使われ, yの各要素がおおむね次式を満たすようにする (i = 0 〜 n-1).   abs(局所誤差) <= rtol*abs(y[i]) + atol atol = 0 と設定すると純粋に相対誤差テストとなる. 誤差が約1.0e-8以下の相対精度を必要とするならば, 通常はderkfを使わないほうがよい. DEPACのdeabmルーチンは厳しい精度を必要とする場合により効率的である. [in] mode (省略可) 動作モード. (省略時 = 0) = 0: toutで戻る (インターバル・モード). = 1: ステップごとに戻る (中間結果出力モード). = 2: ステップごとに戻る (中間結果出力モード). 密出力を可能にする. すなわち, info1 = 2 (最後のステップでは info1 = 1) で戻ったときに derkf_int ルーチンを使って直近のステップの区間内で補間値を求めることができる. (他の値であれば mode = 0 とみなす.) 出典 SLATEC (DEPAC) 参考文献 W H Enright et al. "Interpolants for Runge-Kutta Formulas" ACM Transactions on Mathematical Software Vol.12, No.3, 1986, pp.193-218 使用例 (1) 次の常微分方程式の初期値問題を解く. dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) def f(n, t, y, yp): yp[0] = -2*y[0] + y[1] - cos(t) yp[1] = 2*y[0] - 3*y[1] + 3*cos(t) - sin(t) def TestDerkf(): n = 2 wsave = np.empty(7*n + 20) iwsave = np.empty(20, dtype=int) t = 0 y = np.array([1.0, 2.0]) tfinal = 10.0 tprint = 1.0 info = 0 while True: tout = t + tprint t, info = derkf(info, n, f, t, y, tout, wsave, iwsave) print(t, y) if t >= tfinal or info != 1: break if info != 1: print('Error: info =', info) 実行結果 >>> TestDerkf() 1.0 [0.36787944 0.90818175] 2.0 [ 0.13533528 -0.28081155] 3.0 [ 0.04978707 -0.94020543] 4.0 [ 0.01831564 -0.63532798] 5.0 [0.00673795 0.29040013] 6.0 [0.00247875 0.96264904] 7.0 [0.00091188 0.75481414] 8.0 [ 0.00033546 -0.14516457] 9.0 [ 1.23409832e-04 -9.11006852e-01] 10.0 [ 4.53999600e-05 -8.39026129e-01] 使用例 (2) 次の常微分方程式の初期値問題を解く(密出力). dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) def f(n, t, y, yp): yp[0] = -2*y[0] + y[1] - cos(t) yp[1] = 2*y[0] - 3*y[1] + 3*cos(t) - sin(t) def TestDerkf_2(): n = 2 wsave = np.empty(9*n + 20) iwsave = np.empty(20, dtype=int) y_temp = np.empty(n) t = 0 y = np.array([1.0, 2.0]) tfinal = 10.0 tprint = 1.0 tout = tprint info = 0 while True: t, info = derkf(info, n, f, t, y, tfinal, wsave, iwsave, mode = 2) if info != 1 and info != 2: print('Error: info =', info) break while t >= tout: derkf_int(n, tout, y_temp, wsave) print(tout, y_temp) tout = tout + tprint if t >= tfinal: break 実行結果 >>> TestDerkf_2() 1.0 [0.36787944 0.90818175] 2.0 [ 0.13533528 -0.28081155] 3.0 [ 0.04978707 -0.94020543] 4.0 [ 0.01831564 -0.63532798] 5.0 [0.00673795 0.29040013] 6.0 [0.00247875 0.96264904] 7.0 [0.00091188 0.75481414] 8.0 [ 0.00033546 -0.14516457] 9.0 [ 1.23409855e-04 -9.11006852e-01] 10.0 [ 4.53999601e-05 -8.39026129e-01]