derkfa dopn43

derkfa() def derkfa ( info  , n  , f  , t  , y  , tout  , tend  , wsave  , iwsave  , mode  = 2, rtol  = 1.0e-10, atol  = 1.0e-10  ) 常微分方程式の初期値問題 (5(4)次 ルンゲ・クッタ・フェールベルグ法) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本プログラムは, 5(4)次 ルンゲ・クッタ・フェールベルグ法プログラム DDERKF (文献 (1)) に Enrightら (文献 (2)) の補間式を使用した密出力機能を追加して書き直したものである. 戻り値 (t1, info1) t1 (float): 積分終了時には t = tend, 途中結果を返すために戻ったときには mode の設定により t = tout または t = (直前のステップの終点の値) を返す. info1 (int): リターンコード. = 0: 正常終了. tend までの積分が完了した. < 0: (-info)番目の入力パラメータの誤り. = 1: mode = 1 の中間結果出力のために戻った. 再呼び出しすることより次のステップまで進むことができる. = 2: mode = 2 の中間結果出力のために t = tout において戻った. tout を再設定して再呼び出しすることができる. = 3: mode = 3 の中間結果出力のために t = tout において戻った. tout を再設定して再呼び出しすることができる. = 11: (エラー) 最大ステップ数を超えた. = 12: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. = 13: (エラー) スティフな問題の可能性がある. = 14: (エラー) 誤差テストに失敗した. 誤差許容値が厳しすぎる可能性がある. mode = 1, 2 または 3 における再呼び出し時には t = t1, info = info1 とすること. 引数 [in] info 制御コード. = 0: 最初の呼び出し時(新たに問題を開始する場合)には info = 0 と設定すること. 全ての変数の初期化を行ってから計算を開始する. = 1, 2, 3: info1 = 1, 2 または 3 で戻った場合, 計算を継続するために tout を変更し info = info1 と設定して再呼び出しすることができる. [in] n 微分方程式の数. (n >= 1) [in] f 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. def f(n, t, y, yp): dy/dt を求め yp[i] に入れる (i = 0 〜 n - 1). ただし, n は方程式の数, dy/dt (= f(t, y)) は与えられた t および y における微分の計算値である. [in] t 独立変数 t の現在の値 (1回目の呼び出しにおいては t の初期値 = t0). 本プログラムは t から tend まで積分を行う. mode の設定によって t1 における途中結果を返すために tout においてまたはステップごとに戻ることができる. [in,out] y Numpy ndarray (1次元配列, float, 長さ n) 従属変数 y を表す. [in] t の現在の値における y の値 (1回目の呼び出しにおいては t0 における y の値 = y0). [out] t1 における y の値 (数値解). [in] tout mode = 2 または 3 の場合に, 途中結果の確認／出力を行う t を表す. mode = 0 または 1 では tout は参照されない. mode = 2 または 3 では t = tout における y を求め, それぞれ info1 = 2 または 3 として戻る. 続いて新たな tout における解を求めるために積分を継続したい場合には, info = info1, tout = 新たな tout と設定し, 他の変数を変更せずに再呼び出しすることができる. t < tout <= tend でなければならない (ただし, 後退方向に積分を行う場合には tend <= tout < t でなければならない). 先に tend に達した場合にはその時点で積分を終了し t1 = tend, info1 = 0 として戻る. [out] wsave Numpy ndarray (1次元配列, float, 長さ 9*n + 40) データ領域 (info = 0 として呼び出してから一連の計算の間変更してはならない). [out] iwsave Numpy ndarray (1次元配列, int, 長さ 40) 整数データ領域 (info = 0 として呼び出してから一連の計算の間変更してはならない). [in] mode (省略可) 動作モード. (省略時 = 2) 本プログラムは t から tend までの積分を行うが, 中間結果の確認／出力の方法により4つの動作モードが提供される. どのモードでも t = tend に達したときには積分を終了し, info1 = 0 を返す. = 0: tend まで戻らない. tout は無視される. = 1: 成功したステップごとに戻る (info1 = 1 を返す). tout は無視される. t1 にはそのステップの終了点を返す. 最終ステップでは t1 = tend となるようにステップ幅が調整される. = 2: 積分途中で tout における中間結果を返すために戻る (t1 = tout, info1 = 2 を返す). 続けて次の tout における解を求めるために積分を継続するためは, tout を新たな値に変えて再度呼び出しを行うことができる. tout 直前のステップでは t1 = tout となるようにステップ幅が調節される. = 3: 積分途中で tout における中間結果を返すために戻る (t1 = tout, info1 = 3 を返す). 続けて次の tout における解を求めるために積分を継続するためは, tout を新たな値に変えて再度呼び出しを行うことができる. mode = 2 と異なり tout における y の値は補間により求められ, tout 直前のステップでのステップ幅の調節は行わず mode = 1 と同じステップを踏む. [in] rtol (省略可) 求める解の精度を指定する相対誤差許容値. (省略時 = 1.0e-10) 許容値は atol と共に各ステップにおける局所誤差テストに用いられ, y[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. abs(局所誤差) <= rtol*abs(y[i]) + atol. rtol と atol が同時に 0 であってはならない. [in] atol (省略可) 求める解の精度を指定する絶対誤差許容値. (省略時 = 1.0e-10) 許容値は rtol と共に各ステップにおける局所誤差テストに用いらる (上記 rtol 参照). 文献 (1) SLATEC (DEPAC) (2) W. H. Enright, et al.: "Interpolants for Runge-Kutta Formulas", ACM Transactions on Mathematocal Software, Vol.12, No.3. (1986) 使用例 次の常微分方程式の初期値問題を解く. 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 TestDerkfa(): n = 2 wsave = np.empty(9*n + 40) iwsave = np.empty(40, dtype=int) t = 0 y = np.array([1.0, 2.0]) tend = 10.0 info = 0 while True: tout = t + 1 t, info = derkfa(info, n, f, t, y, tout, tend, wsave, iwsave) print(t, y) if t >= tend or info < 0 or info > 10: break if info != 0: 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]