derkfa dopn43

dopn43() def dopn43 ( info  , n  , f2  , t  , y  , yp  , tout  , tend  , wsave  , iwsave  , mode  = 2, rtol  = 1.0e-10, atol  = 1.0e-10  ) 常微分方程式の初期値問題 (4(3)次ルンゲ・クッタ・ニュストレム法) (2階微分方程式用) 目的 本プログラムは2階常微分方程式の初期値問題 d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y0' の解を求める. ただし, y および y' は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0, y0 および y'0 はそれぞれ t, y および y' の既知の初期値である. 本プログラムは 4(3)次ルンゲ・クッタ・ニュストレム法プログラムである. 係数 (RKN4(3)4FM) は文献(1)を, 密出力のための補間式は文献(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: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. 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] f2 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. def f2(n, t, y, ypp): d2y/dt2 を求め ypp[i] に入れる (i = 0 〜 n - 1). ただし, n は方程式の数, d2y/dt2 (= f2(t, y)) は与えられた t および y における2次微分の計算値である. [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,out] yp Numpy ndarray (1次元配列, float, 長さ n) [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, 長さ 8*n + 40) データ領域 (info = 0 として呼び出してから一連の計算の間変更してはならない). [out] iwsave Numpy ndarray (1次元配列, int, 長さ 30) 整数データ領域 (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 (省略可) 求める解の精度を指定する相対誤差許容値. (rtol >= 0) (省略時 = 1.0e-10) 許容値は atol と共に各ステップにおける局所誤差テストに用いられ, y[] の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. abs(局所誤差) <= rtol*abs(y[i]) + atol および rtol*abs(yp[i]) + atol. rtol と atol が同時に 0 であってはならない. [in] atol (省略可) 求める解の精度を指定する絶対誤差許容値. (atol >= 0) (省略時 = 1.0e-10) 許容値は rtol と共に各ステップにおける局所誤差テストに用いらる (上記 rtol 参照). 文献 (1) J. R. Dormand, M. E. A. El-Mikkawy, P. J. Prince, "Families of Runge-Kutta-Nystrom Formulae", IMA J. of Numerical Analysis, 7, 235-250 (1987). (2) J. R. Dormand, P. J. Prince "Runge-Kutta-Nystrom Triples", Comput. Math. Applic. Vol. 13, No. 12, 937-949 (1987). 使用例 次の2階常微分方程式の初期値問題を解く. y1'' = -y1/R y2'' = -y2/R ただし, R = (y1^2 + y2^2)^(3/2)/(π/4)^2 (t = 0 において y1 = 0.75, y2 = 0, y1' = 0, y2' = (π/4)√(1.25/0.75)) def f(n, t, y, ypp): alfa = pi/4 r = (y[0]**2 + y[1]**2)**(3/2)/alfa**2 ypp[0] = -y[0]/r ypp[1] = -y[1]/r def TestDopn43(): n = 2 wsave = np.empty(8*n + 40) iwsave = np.empty(30, dtype=int) alfa = pi/4 ecc = 0.25 t = 0.0 y = np.array([1 - ecc, 0]) yp = np.array([0, alfa*sqrt((1 + ecc)/(1 - ecc))]) tend = 12.0 info = 0 while True: tout = t + 1 t, info = dopn43(info, n, f, t, y, yp, tout, tend, wsave, iwsave) print(t, y, yp) if t >= tend or info < 0 or info > 10: break if info != 0: print('Error: info =', info) 実行結果 >>> TestDopn43() 1.0 [0.29441755 0.81217852] [-0.7625959 0.47923261] 2.0 [-0.49029978 0.939875 ] [-0.71918028 -0.17238216] 3.0 [-1.05403152 0.5757061 ] [-0.38882952 -0.50909957] 4.0 [-1.25000003e+00 2.61948925e-08] [-3.95363884e-08 -6.08366802e-01] 5.0 [-1.0540316 -0.57570606] [ 0.38882943 -0.5090996 ] 6.0 [-0.49029994 -0.93987504] [ 0.7191802 -0.17238227] 7.0 [ 0.29441735 -0.81217875] [0.76259596 0.47923233] 8.0 [ 7.50000100e-01 -4.82957596e-07] [5.23304832e-07 1.01394459e+00] 9.0 [0.2944181 0.81217828] [-0.76259567 0.47923308] 10.0 [-0.49029923 0.93987522] [-0.71918045 -0.17238178] 11.0 [-1.05403122 0.5757066 ] [-0.38882984 -0.50909938] 12.0 [-1.25000009e+00 6.15458176e-07] [-4.20956957e-07 -6.08366803e-01]