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

Odex2() Sub Odex2 ( N As  Long, F2 As  LongPtr, T As  Double, Y() As  Double, Yp() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Solout2 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 Km As  Long = 0, Optional Nsequ As  Long = 0, Optional Iderr As  Long = 0, Optional Mudif As  Long = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Fac3 As  Double = 0, Optional Fac4 As  Double = 0, Optional Safe1 As  Double = 0, Optional Safe2 As  Double = 0, Optional Safe3 As  Double = 0, Optional Cnt As  Long = 0  ) 常微分方程式の初期値問題 (補外法) (2階微分方程式用) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは2階の常微分方程式の初期値問題 d2y/dt2 = f(t, y), ただし t = t0 において y = y0, y' = y'0 の解を求める. ただし, t0, y0およびy'0は既知でそれぞれt, yおよびy'の初期値である. 上の方程式が連立微分方程式であれば, yおよびy'はベクトルで表される. Odex2はシュテルマー公式にもとづく補外法のプログラムである. ステップ制御, 次数選択および密出力機能を持っている. 詳細については下記参考文献を参照せよ. 引数 [in] N 微分方程式の数. (N >= 1) [in] F2 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F2(N As Long, T As Double, Y() As Double, Yp2() As Double) Yp2(i) = TおよびY()における2階微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Yp2()は与えられたTおよびY()における2階微分係数の計算値, すなわち, Yp2(i) = d2yi/dt2 = fi(T, Y(0), ..., Y(N-1)) (i = 0 〜 N-1). Yp2()以外の変数を変更してはならない. [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数 t の初期値. [out] ステップ終了時の最終値(通常Toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなToutで Info = 1 として再呼び出しすることにより, 新たな点まで積分を継続することができる. [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 = 2N)のどちらでもよい. LRTol = 2, ..., 2N-1の場合, LRTol = 1 とみなす. LRTol > 2N の場合, LRTol = 2N とみなす. LRTol = 2N であっても LATol = 1 であれば LRTol = 1 とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()およびYp()の各要素がおおむね次式を満たすようにする.   スカラーの場合：     abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)   ベクトルの場合 (i = 0 〜 N-1)：     abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. [in] ATol() 配列 ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)またはベクトル(LATol = 2N)のどちらでもよい. LATol = 2, ..., 2N-1の場合, LATol = 1 とみなす. LATol > 2N の場合, LATol = 2N とみなす. LATol = 2N であっても LRTol = 1 であれば LATol = 1 とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()およびYp()の各要素がおおむね次式を満たすようにする.   スカラーの場合：     abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N)   ベクトルの場合 (i = 0 〜 N-1)：     abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i)     abs(Yp(i)の局所誤差) <= RTol(i+N)*abs(Y(i)) + ATol(i+N) 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: 正常終了. = 2: Solout2による中断 (正常終了). = 11: 最大ステップ数を超えた. [in] Solout2 (省略可) 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. (省略時 = NullPtr) Sub Solout2(Nr As Long, Told As Double, T As Double, Y() As Double, Yp() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long) Nr番目のステップTでのY()およびYp()の値が渡されるのでこれを出力する. Toldは前回のTの値, Nは方程式の次数である. Irtrnの値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2 である. Irtrnを使用して積分を中断することができる. Solout2でIrtrnを負の値に設定して戻ると積分を中断して Info = 2 で終了する. Solout2の中で解を変更した場合には Irtrn = 3 と設定して戻ること. 制御情報 RCont および ICont を使用して密出力を行うことができる. 区間[Told, T]内の任意の点T2での解Y(i) (0 <= i <= N-1)を次の関数呼び出しにより得ることができる. Y(i) = Contx2(i, T2, RCont, ICont); End Sub これを省略した場合(Solout2 = NullPtrの場合), 中間結果出力を行わない. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] Hinit (省略可) ステップ幅の初期値. (省略時 = 0.0001) H = 1/||f'|| であり, 通常 0.1〜0.001 が適切とされる. (|Hinit| < 0.0001 であれば Hinit = 0.0001 とみなす) [in] Hmax (省略可) ステップ幅の最大値. (省略時 = Tout - T) (Hmax = 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 許されるステップ数の最大値. (省略時 = 100000) (MaxIter <= 0 であれば省略時の既定値とみなす) [in] Km (省略可) 補外表の最大列数. (Km >= 3) (省略時 = 9) (Km < 3 であれば省略時の既定値とみなす) [in] Nsequ (省略可) ステップ数列の選択. (省略時 = 1 (Solout2不使用時), 4 (Solout2使用時)) = 1: 2, 4, 6, 8, 10, 12, 14, 16, ... = 2: 2, 4, 8, 12, 16, 20, 24, 28, ... = 3: 2, 4, 6, 8, 12, 16, 24, 32, ... = 4: 2, 6, 10, 14, 18, 22, 26, 30, ... Solout2使用時には1〜3は指定できない. (上記以外の値, あるいは, Solout2使用時に1〜3であれば省略時の既定値とみなす) [in] Iderr (省略可) 密出力の式における誤差推定. (省略時 = 0) = 0: 行う = 1: 行わない Solout2を使用しない場合, Iderrの指定にかかわらず誤差推定は行わない. (上記以外の値であれば Iderr = 1 とみなす) [in] Mudif (省略可) 補外式の次数を決めるパラメータ. (1 <= Mudif <= 8) (省略時 = 6) (Mudif < 1 または Mudif > 8 であれば省略時の既定値とみなす) [in] Fac1 (省略可) [in] Fac2 (省略可) ステップ幅選択パラメータ. (省略時: Fac1 = 0.02, Fac2 = 4) j番目の要素が Facmin/Fac2 <= j番目のHnew/Hold <= 1/Facmin となるように新しいステップ幅が選ばれる. ただし, Facmin = Fac1^(1/(2*j - 1)) である. (Fac1 = 0 あるいは Fac2 = 0 であればそれぞれ省略時の既定値とみなす) [in] Fac3 (省略可) [in] Fac4 (省略可) 次数選択パラメータ. (省略時: Fac3 = 0.8, Fac4 = 0.9)   W(k-1) <= W(k)*Fac3 : 次数を減らす   W(k) <= W(k-1)*Fac4 : 次数を増やす W(k) = A(k)/H(k)は単位ステップあたりの仕事量である. (A(k)は関数値計算回数, H(k)はステップ幅を表す) (Fac3 = 0 あるいは Fac4 = 0 であればそれぞれ省略時の既定値とみなす) [in] Safe1 (省略可) [in] Safe2 (省略可) ステップ幅制御の安全係数. (省略時: Safe1 = 0.65, Safe2 = 0.94) (Safe1 = 0 あるいは Safe2 = 0 であればそれぞれ省略時の既定値とみなす) [in] Safe3 (省略可) ステップ幅を減らす際の係数. (省略時 = 0.5) (Safe3 <= 0 であれば省略時の既定値とみなす. Safe3 > 1 であれば 1 とみなす) [in] Cnt (省略可) Neval, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0) = 0: 本ルーチンが呼び出されるたびにリセットする. <> 0: Info = 0 で呼び出されたときにだけリセットする. 出典 E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993) 邦訳: 「常微分方程式の数値解法Ⅰ 基礎編」スプリンガージャパン (2007) 使用例 (1) 次の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_Odex2() 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 Odex2(N, AddressOf F2, T, Y(), Yp(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Odex2: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 -0.427967245621361 0.863775701003976 2 -1.20572535244133 0.613566455365295 3 -1.4955436794428 8.16675371786414E-02 4 -1.33475968929777 -0.47684609244579 5 -0.700827262135424 -0.848381581815096 6 0.357480600846299 -0.44558418320875 7 -0.118067376945689 0.800372165255964 8 -1.03762003516232 0.730221556218178 9 -1.45981364707677 0.243039748770668 10 -1.42617024854319 -0.326583069133528 使用例 (2) 次の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_Odex2_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tend 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 Call Odex2(N, AddressOf F2, T, Y(), Yp(), Tend, RTol(), ATol(), Info, AddressOf SoloutX2) If Info <> 1 Then Debug.Print "Error in Odex2: Info =", Info End Sub Sub SoloutX2(Nr As Long, Told As Double, T As Double, Y() As Double, Yp() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long) Dim Y0 As Double, Y1 As Double Static Tout As Double If Nr = 1 Then Tout = 1 While T >= Tout Y0 = Contx2(0, Tout, RCont, ICont) Y1 = Contx2(1, Tout, RCont, ICont) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End Sub 実行結果 1 -0.427967245557026 0.863775700969804 2 -1.20572535233059 0.613566455161914 3 -1.49554367919972 8.16675368902836E-02 4 -1.33475968872957 -0.476846092763801 5 -0.700827261091603 -0.848381581795345 6 0.357480601918728 -0.445584181833266 7 -0.118067378360736 0.800372165912016 8 -1.03762003630738 0.730221557185247 9 -1.45981364907474 0.243039751005974 10 -1.42617025292365 -0.326583065997367