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

Radau5_r() Sub Radau5_r ( N 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 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 MaxNit As  Long, Optional StartN As  Long, Optional Nind1 As  Long, Optional Nind2 As  Long, Optional Nind3 As  Long, Optional Pred As  Long, Optional Hess As  Long, Optional Safe As  Double, Optional Thet As  Double, Optional FNewt As  Double, Optional Quot1 As  Double, Optional Quot2 As  Double, Optional Facl As  Double, Optional Facr As  Double, Optional Cnt As  Long  ) 常微分方程式の初期値問題 (5次の陰的ルンゲ・クッタ法 (ラダウIIA法)) (リバースコミュニケーション版) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは1階のスティフな常微分方程式 (あるいは微分代数方程式)の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Mは質量マトリックスである. 方程式は線形陰的(M = I(単位行列) でない)でも, 陽的(M = I)でもよい. Radau5は5次のラダウIIA法に基づく陰的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. Radau5_rはRadau5のリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [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) = -3: パラメータ Y() の誤り. = -5: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -6: パラメータ ATol() の誤り. (ATol(i) < 0) = -7: パラメータ Cont() の誤り. = -8: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = -10: パラメータ YY() の誤り. = -11: パラメータ YYp() の誤り. = -12: パラメータ YYpd() の誤り. = -14: パラメータ IRev の誤り. (IRev <> 0, 1, 2, 4, 5) = -20: パラメータ Mljac あるいは Mlmas の誤り. (Mlmas > Mljac) = -21: パラメータ Mujac あるいは Mumas の誤り. (Mumas > Mujac) = -27: パラメータ M1 の誤り. (M1 < 0) = -28: パラメータ M1 または M2 の誤り. (M2 < 0 または M1 + M2 > N) = -34: パラメータ Nind1, Nind2 あるいは Nind3 の誤り. (Nind1 + Nind2 + Nind3 <> N) = -41: パラメータ FNewt の誤り. (FNewt が小さすぎる) = 1: 正常終了. = 2: Irtrn による中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: 行列が特異になった. [out] TT IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のTの値を返す. [out] YY() 配列 YY(LYY - 1) (LYY >= N) IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のYの値を返す. [in] YYp() 配列 YYp(LYYp - 1) (LYYp >= N) IRev = 1: T(= TT)およびY(= YY())における微分係数, すなわち, YYp(i) = dyi/dt = fi(TT, YY(0), ..., YY(N-1)) (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()以外の変数を変更してはならない. = 4: 質量マトリックスMをYYpd()に設定する. YYpd()以外の変数を変更してはならない. = 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). [in] Iout (省略可) 計算の途中結果を出力するかどうか指定する. (省略時 = 0) = 0: 出力しない. (IRev = 5で戻らない) = 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でないRadau5のSoloutの呼び出しに相当する. 対応する情報は次のとおりである. Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY(). Cont()はSoloutと同様に密出力のために使用する. Y(i) = Contr5_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] 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] MaxNit (省略可) 陰形式の場合の各ステップにおける解の計算のためのニュートン法の最大反復回数. (省略時 = 7) (MaxNit <= 0 であれば省略時の既定値とみなす) [in] StartN (省略可) ニュートン法の出発値. (省略時 = 0) = 0: 外挿した解を使う. <> 0: 0 を使う. ニュートン法の収束に問題がある場合には後者が推奨される. [in] Nind1,Nind2,Nind3 (省略可) 指数1, 指数2, 指数3変数それぞれの次数. (Nind1 > 0, Nind1 + Nind2 + Nind3 = N) (省略時: Nind1 = N, Nind2 = 0, Nind3 = 0) 指数 > 1 の微分代数方程式(DAE)に関するパラメータである. 関数を定義するサブルーチンは, 指数1, 2および3の変数がこの順番に現れるように作成されていなけらばならない. 常微分方程式(ODE)の場合, Nind1は方程式の次数に一致する. [in] Pred (省略可) ステップ幅戦略の選択. (省略時 = 1) = 1: モデル予期コントローラ(グスタフソン). = 2: 伝統的なステップ幅コントロール. (上記以外の値であれば省略時の既定値とみなす) [in] Hess (省略可) ヤコビ行列をヘッセンベルグ形に変換するかどうか指定する. (省略時 = 0) = 0: 変換しない. (本機能はサポートされない. 0 以外の値であれば 0 を指定したものとみなす.) [in] Safe (省略可) ステップ幅推定時の安全係数. (0.001 < Safe < 1) (省略時 = 0.9) (Safe <= 0.001 あるいは Safe >= 1 であれば省略時の既定値とみなす) [in] Thet (省略可) ヤコビ行列を再計算するかどうか指定する. (Thet < 1) (省略時 = 0.001) ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい(例えば 0.1). 小規模な場合には小さくすればよい(例えば 0.001). 負の値にすると, ステップが採用されるたびに強制的にヤコビ行列が再計算される. (Thet = 0 あるいは Thet >= 1 であれば省略時の既定値とみなす) [in] FNewt (省略可) ニュートン法の停止基準. (通常 FNewt < 1) (省略時 = RTol(0)より適切な値を設定) 小さな値にすると遅くなるがより安全になる. (FNewt = 0 であれば省略時の既定値とみなす) [in] Quot1,Quot2 (省略可) Quot1 < Hnew/Hold < Quot2 であればステップ幅は変更されない. (Quot1 <= 1, Quot2 >= 1) (省略時: Quot1 = 1, Quot2 = 1.2) 小規模な場合 Quot1 = 1, Quot2 = 1.2, 大規模な場合 Quot1 = 0.99, Quot2 = 2 とすればよい. (Quot1 = 0 または Quot1 > 1, Quot2 = 0 または Quot2 < 1 であればそれぞれ省略時の既定値とみなす) [in] Facl,Facr (省略可) ステップ幅選択パラメータ. (Facl <= 1, Facr >= 1) (省略時: Facl = 0.2, Facr = 8) Facl < Hnew/Hold < Facr となるようにステップ幅が選ばれる. (Facl = 0 または Facl > 1, Facr = 0 または Facr < 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_Radau5_r() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Iout As Long, 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, Irtrn As Long, IRev As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Radau5_r(N, 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 Radau5_r: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441187316 0.90818174710239 2 0.135335283255029 -0.280811553607639 3 4.97870683911592E-02 -0.940205428513689 4 1.83156389200067E-02 -0.635327984370112 5 6.73794702805148E-03 0.290400132651525 6 2.47875220047626E-03 0.96264903967345 7 9.11881982753585E-04 0.754814137453754 8 3.35462642879546E-04 -0.145164572401227 9 1.23409800092344E-04 -0.911006854141312 10 4.53999251935239E-05 -0.839026134868127 使用例 (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_Radau5_r_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Iout As Long, 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, Irtrn As Long, IRev As Long, Y0 As Double, Y1 As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Tout = 1 Info = 0 IRev = 0 Do Call Radau5_r(N, 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 = Contr5_r(0, Tout, Cont()) Y1 = Contr5_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 Radau5_r: Info =", Info End Sub 実行結果 1 0.36787944514433 0.90818175052153 2 0.135335281095574 -0.280811550970382 3 4.97870683819513E-02 -0.940205439380229 4 1.83156414146801E-02 -0.635327955429466 5 6.73794660800666E-03 0.29040008345163 6 2.47875304716379E-03 0.962649009570966 7 9.11881374205137E-04 0.754814090253208 8 3.35462420859941E-04 -0.145164626236664 9 1.23409926575974E-04 -0.911006958450603 10 4.53999332681047E-05 -0.839026159221514