Contex Contex_r Contr5 Contr5_r Contra Contra_r Contro Contro_r Contrp Contrp_r Dassl Dassl_r Radau Radau5 Radau5_r Radau_r Radaup Radaup_r Rodas Rodas_r Seulex Seulex_r

Rodas_r() Sub Rodas_r ( N As  Long, Ifcn As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Cont() As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, YYpd() As  Double, Irtrn As  Long, IRev As  Long, Optional Iout As  Long, Optional Ijac As  Long, Optional Mljac As  Long = -1, Optional Mujac As  Long, Optional Idfx As  Long, Optional Imas As  Long, Optional Mlmas As  Long = -1, Optional Mumas As  Long, Optional Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional M1 As  Long, Optional M2 As  Long, Optional Hinit As  Double, Optional Hmax As  Double, Optional MaxIter As  Long, Optional Meth As  Long, Optional Pred As  Long, Optional Safe As  Double, Optional Fac1 As  Double, Optional Fac2 As  Double, Optional Cnt As  Long  ) 常微分方程式の初期値問題 (4(3)次 ローゼンブロック法) (リバースコミュニケーション版) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Mは質量マトリックスである. 方程式は線形陰的(M = I(単位行列) でない)でも, 陽的(M = I)でもよい. Rodasは4(3)次のローゼンブロック法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. Rodas_rはRodasのリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [in] Ifcn f()がtに依存するかどうかを示す. = 0: tに依存しない dy/dt = f(y) (自励系). = 1: tに依存する dy/dt = f(y, t) (非自励系). [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] 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になってはいけない. [out] Cont() 配列 Cont(LCont - 1) (LCont >= 4*N) 密出力のための制御情報. (Iout = 0 のときには参照されない) [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: パラメータ N の誤り. (N < 1) = -4: パラメータ Y() の誤り. = -6: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -7: パラメータ ATol() の誤り. (ATol(i) < 0) = -8: パラメータ Cont() の誤り. = -9: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = -11: パラメータ YY() の誤り. = -12: パラメータ YYp() の誤り. = -13: パラメータ YYpd() の誤り. = -15: パラメータ IRev の誤り. (IRev <> 0, 1, 2, 3, 4, 5) = -22: パラメータ Mljac あるいは Mlmas の誤り. (Mlmas > Mljac) = -23: パラメータ Mujac あるいは Mumas の誤り. (Mumas > Mujac) = -29: パラメータ M1 の誤り. (M1 < 0) = -30: パラメータ M1 または M2 の誤り. (M2 < 0 または M1 + M2 > N) = 1: 正常終了. = 2: Irtrn による中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: 行列が特異になった. [out] TT IRev = 1, 2, 3: 微分係数あるいはヤコビ行列を求めてYYp()に入れるべき点の独立変数tの値. [out] YY() 配列 YY(LYY - 1) (LYY >= N) IRev = 1, 2, 3: 微分係数あるいはヤコビ行列を求めてYYp()に入れるべき点の従属変数yの値. [in] YYp() 配列 YYp(LYYp - 1) (LYYp >= N) IRev = 1: TTおよびYY()における微分係数, dyi/dt = fi(TT, YY()) (i = 0 〜 N-1) を再呼び出し時に入力する. IRev = 3: TTおよびYY()におけるfのtについての偏微分, dfi(TT, YY())/dt (i = 0 〜 N-1) を再呼び出し時に入力する. [in] YYpd() 配列 YYpd(LYYpd1 - 1, LYYpd2 - 1) (LYYpd1 >= max(Ljac, Lmas), LYYpd2 >= N) IRev = 2: TTおよびYY()におけるヤコビ行列(dfi/dyj)を再呼び出し時に入力する.   注 - Mljac = N の場合は通常の2次元配列(Ljac = N), Mljac < N の場合は帯行列形式(Ljac = Mljac + Mujac + 1)で格納する. Ijac = 0 ならば Ljac = 0 とする. IRev = 4: 質量マトリクスMを再呼び出し時に入力する.   注 - Mlmas = N の場合は通常の2次元配列(Lmas = N), Mlmas < N の場合は帯行列形式(Lmas = Mlmas + Mumas + 1)で格納する. Imas = 0 ならば Lmas = 0 とする. [in,out] Irtrn [in] IRev = 5: 処理を中断したい場合以外変更してはならない. Irtrnを負の値に設定すると処理を中断して Info = 2 で終了する. [out] IRev = 5: 初回, 中間, 最後にIRev = 5 で戻るとき, それぞれ 0, 1, 2 を返す. [in,out] IRev リバースコミュニケーションの制御変数 [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: TTおよびYY()における微分係数の計算しYYp()に設定する. YYp()以外の変数を変更してはならない. = 2: TTおよびYY()におけるヤコビ行列を計算しYYpd()に設定する. YYpd()以外の変数を変更してはならない. = 3: TTおよびYY()におけるfiのtによる偏微分(dfi/dt)を計算しYYpd()に設定する. YYpd()以外の変数を変更してはならない. = 4: 質量マトリックスMをYYpd()に設定する. YYpd()以外の変数を変更してはならない. = 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). [in] Iout (省略可) 計算の途中結果を出力するかどうか指定する. (省略時 = 0) = 0: 出力しない. (IRev = 5で戻らない) = 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でないRodasのSoloutの呼び出しに相当する. 対応する情報は次のとおりである. Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY(). Cont()はSoloutと同様に密出力のために使用する. Y(i) = Contro_r(i, T2, Cont()) (Ioutが上記以外の値であれば省略時の既定値とみなす) [in] Ijac (省略可) ヤコビ行列の計算方法の選択. (省略時 = 0) = 0: ヤコビ行列を有限差分により求める. (IRev = 2で戻らない) = 1: ヤコビ行列を外部でユーザが計算する. (計算が必要なときにIRev = 2で戻る) (他の値であれば省略時の既定値とみなす) [in] Mljac (省略可) ヤコビ行列の下帯幅. (0 <= Mljac <= N) (省略時 = N) Mljac = N の場合, ヤコビ行列はN×N行列として格納される. Mljac < N の場合, ヤコビ行列は帯行列として格納される. (Mljac < 0 または Mljac > N であれば省略時の既定値とみなす) [in] Mujac (省略可) ヤコビ行列の上帯幅. (0 <= Mujac <= N) (省略時 = 0) Mljac = N の場合, Mujacは無視される. (Mujac < 0 または Mujac > N であれば省略時の既定値とみなす) [in] Idfx (省略可) f(t, y)のtについての偏微分の計算方法の選択. (省略時 = 0) = 0: 有限差分を使って求める. (IRev = 3では戻らない) = 1: 外部でユーザーが計算する. (計算が必要なときにIRev = 3で戻る) (Ifcn = 1 の場合に有効) [in] Imas (省略可) 質量マトリクスMの計算方法の選択. (省略時 = 0) = 0: Mは単位行列とする. (IRev = 4で戻らない) = 1: 質量マトリクスMを外部でユーザーが計算する. (計算が必要なときにIRev = 4で戻る) (他の値であれば省略時の既定値とみなす) [in] Mlmas (省略可) 質量マトリクスMの下帯幅. (0 <= Mljac <= N) (省略時 = N) Mlmas = N の場合, MはN×N行列として格納される. Mlmas < N の場合, Mは帯行列として格納される. (Mlmas < 0 または Mlmas > N であれば省略時の既定値とみなす) [in] Mumas (省略可) 質量マトリクスMの上帯幅. (0 <= Mumas <= N) (省略時 = 0) Mlmas = N の場合, Mumasは無視される. (Mumas < 0 または Mumas > N であれば省略時の既定値とみなす) [out] Neval (省略可) 関数評価回数. (ヤコビ行列の計算のためのものは含まない) [out] Njac (省略可) ヤコビ行列評価回数. (有限差分の場合を含む) [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] M1,M2 (省略可) 最初のM1本の方程式が次の形式の場合を考える.   y'(i) = y(i + M2) (i = 0 〜 M1-1) M1はM2の整数倍とし, かつ残りの方程式がy'(M1), ..., y'(N-1)に陽に依存しない場合, パラメータM1とM2を0以外に設定すると効率よく計算される. (M1 > 0, M2 > 0, M1 + M2 <= N) (省略時は M1 = M2 = 0) パラメータを0以外に設定としたときには, ヤコビ行列の非ゼロ要素(行 M1〜N-1)だけを(N - M1)×Nの配列に入れなければならない. また, 質量マトリックスMの(N - M1)次元右下ブロックだけを(N - M1)×(N - M1)の配列に入れなければならない. [in] Hinit (省略可) ステップ幅の初期値. (省略時 = 1.0e-6) H = 1/||f'||である. 初期トランジェントのあるスティフな方程式の場合, 通常は1.0e-3〜1.0e-5 が適切とされる. (Hinit = 0 であれば1.0e-6とする) [in] Hmax (省略可) ステップ幅の最大値. (省略時 = Tout - T) (Hmax = 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 許されるステップ数の最大値. (省略時 = 100000) (MaxIter <= 0 であれば省略時の既定値とみなす) [in] Meth (省略可) 係数の選択. (省略時 = 1) = 1: 出典のP452(邦訳 P426)の方法. = 2: 同じ方法で異なる係数を使用. = 3: Gerd Steinebachの係数を使用した方法. (上記以外の値であれば省略時の既定値とみなす) [in] Pred (省略可) ステップ幅戦略の選択. (省略時 = 1) = 1: モデル予期コントローラ(グスタフソン). = 2: 伝統的なステップ幅コントロール. (上記以外の値であれば省略時の既定値とみなす) [in] Safe (省略可) ステップ幅推定時の安全係数. (0.001 < Safe < 1) (省略時 = 0.9) (Safe <= 0.001 あるいは Safe >= 1 であれば省略時の既定値とみなす) [in] Fac1,Fac2 (省略可) ステップ幅選択パラメータ. (Fac1 <= 1, Fac2 >= 1) (省略時: Fac1 = 0.2, Fac2 = 6) Fac1 < Hnew/Hold < Fac2 となるようにステップ幅が選ばれる. (Fac1 = 0 または Fac1 > 1, Fac2 = 0 または Fac2 < 1 であればそれぞれ省略時の既定値とみなす) [in] Cnt (省略可) Neval, Njac, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0) = 0: 本ルーチンが呼び出されるたびにリセットする. <> 0: Info = 0 で呼び出されたときにだけリセットする. 出典 E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations II. Stiff and differential-algebraic Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1996) 邦訳: 「常微分方程式の数値解法Ⅱ 発展編」スプリンガージャパン (2008) 使用例 (1) 次の常微分方程式の初期値問題(スティフな問題)を解く. dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 1998 * Y(0) - 1999 * Y(1) + 1999 * Cos(T) - Sin(T) End Sub Sub Ex_Rodas_r() Const N = 2 Dim Ifcn As Long, T As Double, Y(N - 1) As Double, Tend As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim Cont(4 * N) As Double Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double Dim Irtrn As Long, IRev As Long, Iout As Long, Tout As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Ifcn = 1 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Rodas_r(N, Ifcn, T, Y(), Tout, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Rodas_r: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441171443 0.908181747041726 2 0.135335283236614 -0.280811553312466 3 4.97870683678661E-02 -0.940205428235504 4 1.83156388887351E-02 -0.635327981976607 5 6.73794699908466E-03 0.290400132463827 6 2.47875217666488E-03 0.962649038829728 7 9.11881965553467E-04 0.754814136310732 8 3.35462627903455E-04 -0.145164571182498 9 1.2340980408788E-04 -0.911006852082675 10 4.53999297639715E-05 -0.839026129149385 使用例 (2) 次の常微分方程式の初期値問題(スティフな問題)を解く(密出力を使用). dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 1998*y1 - 1999*y2 + 1999*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) Sub F2(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 1998 * Y(0) - 1999 * Y(1) + 1999 * Cos(T) - Sin(T) End Sub Sub Ex_Rodas_r_2() Const N = 2 Dim Ifcn As Long, T As Double, Y(N - 1) As Double, Tend As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim Cont(4 * N) As Double Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double Dim Irtrn As Long, IRev As Long, Iout As Long, Tout As Double Dim Y0 As Double, Y1 As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 Ifcn = 1 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Tout = 1 Info = 0 IRev = 0 Do Call Rodas_r(N, Ifcn, T, Y(), Tend, RTol(), ATol(), Cont(), Info, TT, YY(), YYp(), YYpd(), Irtrn, IRev, Iout) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) ElseIf IRev = 5 Then While T >= Tout Y0 = Contro_r(0, Tout, Cont()) Y1 = Contro_r(1, Tout, Cont()) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Rodas_r: Info =", Info End Sub 実行結果 1 0.36787944117149 0.908181746948026 2 0.135335283236613 -0.280811553312504 3 4.97870683678624E-02 -0.940205428229929 4 1.83156388886953E-02 -0.635327981896935 5 6.73794699908206E-03 0.290400132468997 6 2.47875217666092E-03 0.962649038837632 7 9.11881965597003E-04 0.754814136223758 8 3.35462627893874E-04 -0.145164571163287 9 1.23409804039225E-04 -0.911006851985428 10 4.53999297635988E-05 -0.839026129148628