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() Sub Rodas ( N As  Long, F As  LongPtr, Ifcn As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Solout As  LongPtr = NullPtr, Optional Jac As  LongPtr = NullPtr, Optional Mljac As  Long = -1, Optional Mujac As  Long, Optional Dfx As  LongPtr = NullPtr, Optional Mas As  LongPtr = NullPtr, 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)次のローゼンブロック法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(i) = TおよびY()における微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Yp()は与えられたTおよびY()における微分係数の計算値, すなわち, Yp(i) = dyi/dt = fi(T, Y(0), ..., Y(N-1)) (i = 0 〜 N-1). Yp()以外の変数を変更してはならない. Mは通常単位行列であるが, 必要であればMasで定義することができる. [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になってはいけない. [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: パラメータ N の誤り. (N < 1) = -5: パラメータ Y() の誤り. = -7: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -8: パラメータ ATol() の誤り. (ATol(i) < 0) = -9: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = -16: パラメータ Mljac あるいは Mlmas の誤り. (Mlmas > Mljac) = -17: パラメータ Mujac あるいは Mumas の誤り. (Mumas > Mujac) = -23: パラメータ M1 の誤り. (M1 < 0) = -24: パラメータ M1 または M2 の誤り. (M2 < 0 または M1 + M2 > N) = 1: 正常終了. = 2: Soloutによる中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: 行列が特異になった. [in] Solout (省略可) 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. (省略時 = NullPtr) Sub Solout(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, Cont As Double, Irtrn As Long) Nr番目のステップTでのY()の値が渡されるのでこれを出力する. Toldは前回のTの値, Nは方程式の次数である. Irtrnの値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2 である. Irtrnを使用して積分を中断することができる. SoloutでIrtrnを負の値に設定して戻ると積分を中断して Info = 2 で終了する. 制御情報 Cont を使用して密出力を行うことができる. 区間[Told, T]内の任意の点T2での解Y(i) (0 <= i <= N-1)を次の関数呼び出しにより得ることができる. Y(i) = Contro(i, T2, Cont); End Sub これを省略した場合(Solout = NullPtrの場合), 中間結果出力を行わない. [in] Jac (省略可) ヤコビ行列 dfi(t, y)/dyj を計算するユーザー定義サブルーチンで, 次のように定義すること. (省略時 = NullPtr) Sub Jac(N As Long, T As Double, Y() As Double, Ypd() As Double) Ypd(i, j) = TおよびY()におけるdfi/dyjの計算値 (i = 0〜N-1, j = 0〜N-1). End Sub Mljac = N の場合Ypd()は通常の2次元配列, 0 <= Mljac < N の場合帯行列形式である. Ypd()以外の変数を変更してはならない. これを省略した場合(Jac = NullPtrの場合), ヤコビ行列は有限差分により求められる. [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] Dfx (省略可) f(t, y)のtについての偏微分を計算するユーザー定義サブルーチンで, 次のように定義すること. Sub Dfx(N As Long, T As Double, Y() As Double, Fx() As Double) Fx(i) = dfi/dt (i = 0〜N-1) End Sub Ifcn = 1 の場合にのみ呼び出される. これを省略した場合(Dfx = NullPtr の場合), 偏微分は有限差分を使って求められる. [in] Mas (省略可) 質量マトリクスMを与えるユーザー定義サブルーチンで, 次のように定義すること. (省略時 = NullPtr) Sub Mas(N As Long, Am() As Double) Am(i, j) = 質量マトリクスMijの値 (i = 0〜N-1, j = 0〜N-1). End Sub Mlmas = N の場合Am()は通常の2次元配列, 0 <= Mlmas < N の場合帯行列形式である. Am()以外の変数を変更してはならない. これを省略した場合(Mas = NullPtrの場合), 質量マトリクスMは単位行列とみなす. [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以外に設定としたときには, Jacにおいてヤコビ行列の非ゼロ要素(行 M1〜N-1)だけを(N - M1)×Nの配列に入れなければならない. また, Masにおいて質量マトリックス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() Const N = 2 Dim Ifcn As Long, T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long, I As Long 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 Call Rodas(N, AddressOf F2, Ifcn, T, Y(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Rodas: 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_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, I As Long 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 Call Rodas(N, AddressOf F2, Ifcn, T, Y(), Tend, RTol(), ATol(), Info, AddressOf SoloutRo) If Info <> 1 Then Debug.Print "Error in Rodas: Info =", Info End Sub Sub SoloutRo(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, Cont As Double, Irtrn As Long) Dim Y0 As Double, Y1 As Double Static Tout As Double If Nr = 1 Then Tout = 1 While T >= Tout Y0 = Contro(0, Tout, Cont) Y1 = Contro(1, Tout, Cont) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend 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