Contd5 Contd5_r Contd8 Contd8_r Contx1 Contx1_r Contx2 Contx2_r Derkf Derkf_r DerkfInt Dop853 Dop853_r Dopri5 Dopri5_r Doprin Doprin_r Dverk Dverk_r DverkInt Odex Odex2 Odex2_r Odex_r Retard Retard_r Ylag Ylag_r

Doprin() Sub Doprin ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Yp() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Solout As  LongPtr = NullPtr, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional MaxIter As  Long = 0, Optional Safe As  Double = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Cnt As  Long = 0  ) 常微分方程式の初期値問題 (7(6)次ルンゲ・クッタ・ニュストレム法) (2階微分方程式用) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは2階の常微分方程式の初期値問題 d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0 の解を求める. ただし, t0, y0およびy'0は既知で, それぞれt, yおよびy'の初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンはドルマンとプリンスによる7(6)次ルンゲ・クッタ・ニュストレム法のプログラムである. ステップ制御アルゴリズムを持っている. 詳細については下記参考文献を参照せよ. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Ypp() As Double) Ypp(i) = TおよびY()における2次微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Ypp()は与えられたTおよびY()における2次微分係数, すなわち, Ypp(i) = d2yi/dt2 = fi(T, Y(0), ..., Y(N-1)) (i = 0 〜 N-1). [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数 t の初期値. [out] ステップ終了時の最終値(通常Toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなToutで再呼び出しすることにより, 新たな点まで積分を継続することができる. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数Y()の初期値. [out] 最終のT(通常Toutに等しい)において求められた解(の近似値). [in,out] Yp() 配列 Yp(LYp - 1) (LYp >= N) [in] Tの初期値における微分値y'1〜y'nの初期値. [out] 最終のT(通常Toutに等しい)において求められた微分値(の近似値). [in] Tout 解を求めたい点を設定する. 積分を行うのはTについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. [in] RTol() 配列 RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)またはベクトル(LRTol = N)のどちらでもよい. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. LRTol = N であっても LATol = 1 であれば LRTol = 1 とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y(i)の各要素がおおむね次式を満たすようにする (i = 0 〜 LRTol-1).   abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. [in] ATol() 配列 ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)またはベクトル(LATol = N)のどちらでもよい. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. LATol = N であっても LRTol = 1 であれば LATol = 1 とみなす. 許容値は各ステップにおけるローカルな誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 LATol-1).   abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: パラメータ N の誤り. (N < 1) = -4: パラメータ Y() の誤り. = -5: パラメータ Yp() の誤り. = -7: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -8: パラメータ ATol() の誤り. (ATol(i) < 0) = -9: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = 1: 正常終了. = 11: 収束しなかった (ステップ数の最大値を超えたなど). [in] Solout (省略可) 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. Sub Solout(Nr As Long, T As Double, Y() As Double, Yp() As Double, N As Long, Irtrn As Long) Nr番目のステップTでのY()およびYp()の値が渡されるのでこれを出力する. Nは方程式の次数である. Irtrnの値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2 である. Irtrnを使用して積分を中断することができる. SoloutでIrtrnを負の値に設定して戻ると積分を中断して Info = 2 で終了する. End Sub これを省略した場合(Solout = NullPtrの場合), 中間結果出力を行わない. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] Hinit (省略可) ステップ幅の初期値. (省略時 = 初期関数値より自動推定) (Hinit = 0 であれば省略時の既定値とみなす) [in] Hmax (省略可) ステップ幅の最大値. (Hmax >= 0) (省略時 = Tout - T) (Hmax = 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 許されるステップ数の最大値. (MaxIter >= 0) (省略時 = 2000) (MaxIter = 0 であれば省略時の既定値とみなす) [in] Safe (省略可) ステップ幅推定時の安全係数. (0.0001 < Safe < 1) (省略時 = 0.9) (Safe <= 0.0001 または Safe >= 1 であれば省略時の既定値とみなす) [in] Fac1 (省略可) [in] Fac2 (省略可) ステップ幅選択パラメータ. (省略時: Fac1 = 0.2, Fac2 = 10) Fac1 <= Hnew/Hold <= Fac2 となるようにステップ幅が選ばれる. (Fac1=0 あるいは Fac2=0 であればそれぞれ省略時の既定値とみなす) [in] Cnt (省略可) Neval, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0) = 0: 本ルーチンが呼び出されるたびにリセットする. <> 0: Info = 0 で呼び出されたときにだけリセットする. 参考文献 E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations I", Springer Series in Computational Mathematics, Springer-Verlag (1987) J. R. Dormand and P. J. Prince, "New Runge-Kutta algorithms for numerical simulation in dynamical astoronomy", Celestial Mechanics 18, 223?232 (1978) 使用例 次の2階常微分方程式の初期値問題を解く. y1'' = -y1/R y2'' = -y2/R ただし, R = (y1^2 + y2^2)^(3/2) (t = 0 において y1 = 0.5, y2 = 0, y1' = 0, y2' = √3) Sub F2(N As Long, T As Double, Y() As Double, Yp2() As Double) Dim R As Double R = (Y(0) ^ 2 + Y(1) ^ 2) ^ (3 / 2) Yp2(0) = -Y(0) / R Yp2(1) = -Y(1) / R End Sub Sub Ex_Doprin() Const N = 2 Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 0.5: Y(1) = 0: Yp(0) = 0: Yp(1) = Sqr(3) Info = 0 Do Tout = T + 1 Call Doprin(N, AddressOf F2, T, Y(), Yp(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Doprin: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 -0.427967245622589 0.863775701069405 2 -1.20572535251092 0.61356645560506 3 -1.49554367984035 8.16675377144076E-02 4 -1.33475969024023 -0.476846091796985 5 -0.700827263990535 -0.848381581520395 6 0.357480599005646 -0.445584185957294 7 -0.118067374412737 0.800372164683605 8 -1.03762003373262 0.730221558340758 9 -1.45981364767163 0.243039752399323 10 -1.42617025207063 -0.326583065015208