Debdf Debdf_r Radaua Radaua_r Rodasa Rodasa_r Seulexa Seulexa_r

Radaua() Sub Radaua ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, RTol() As  Double, ATol() As  Double, Mode As  Long, Info As  Long, Optional Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional Ndec As  Long, Optional Nsol As  Long, Optional Fjac As  LongPtr = NullPtr, Optional Mljac As  Long = -1, Optional Mujac As  Long, Optional Fmas As  LongPtr = NullPtr, Optional Mlmas As  Long = -1, Optional Mumas As  Long, Optional Hes As  Long, Optional MaxIter As  Long, Optional Nit1 As  Long, Optional Startn As  Long, Optional Nind1 As  Long, Optional Nind2 As  Long, Optional Nind3 As  Long, Optional Pred As  Long, Optional M1 As  Long, Optional M2 As  Long, Optional Nsmax As  Long, Optional Nsmin As  Long, Optional Ns As  Long, Optional Cnt As  Long, Optional Hinit As  Double, Optional Hmax As  Double, Optional Thet As  Double, Optional Facl As  Double, Optional Facr As  Double, Optional Safe As  Double, Optional Quot1 As  Double, Optional Quot2 As  Double, Optional Vitu As  Double, Optional Vitd As  Double, Optional Hhou As  Double, Optional Hhod As  Double  ) 常微分方程式の初期値問題 (5, 9, 13次 可変次数陰的ルンゲ・クッタ法 (ラダウIIA法)) 目的 本プログラムは1階のスティフな常微分方程式 (あるいは微分代数方程式) の初期値問題 M * dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. M は質量マトリックスである. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本プログラムは, 5, 9, 13次 可変次数陰的ルンゲ・クッタ法 (ラダウIIA法) プログラム RADAU (文献 (1)) を書き直したものである. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の関数値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double) dy/dt を求め Yp() に入れる. End Sub ただし, N は方程式の数, dy/dt (= f(t, y)) は与えられた T および Y() における微分の計算値である. [in,out] T 独立変数 t を表す. 本プログラムは t の初期値から Tend までの積分を行う. Mode の設定により, 途中結果を返すために Tout またはステップごとに戻ることができる. [in] t の初期値 t0. [out] 積分終了時には T = Tend, 途中結果を返すために戻ったときには Mode の設定により T = Tout または T = 直前のステップの終点の値 を返す. [in,out] Y() 配列 Y(LY - 1) (LY >= N) 従属変数 y を表す. [in] t の初期値 t0 における y の初期値 y0. [out] t = T における y の値 (数値解). [in] Tout Mode = 2 または 3 の場合に, 途中結果の確認／出力を行う t を表す. Mode = 0 または 1 では Tout は参照されない. Mode = 2 または 3 では t = Tout における y を求め, それぞれ Info = 2 または 3 として戻る. 続いて新たな Tout における解を求めるために積分を継続したい場合には, 新たな Tout に変更して(Info を含む)他の変数を変更せずに再度呼び出しを行うことができる. T < Tout <= Tend でなければならない (ただし, 後退方向に積分を行う場合には Tend <= Tout < T でなければならない). 先に Tend に達した場合にはその時点で積分を終了し T = Tend, Info = 0 として戻る. [in] Tend 積分を終了する点を表す. Tend に達すると積分を終了し, T = Tend, Info = 0 として戻る. 積分を行うのは前進方向 (Tend > T) でも後退方向 (Tend < T) でもよい. [in] Rtol() 配列 Rtol(LRtol - 1) (LRtol >= 1) (Rtol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値を表す. 本パラメータは LRtol および LAtol の値によりスカラーまたは配列を選択できる. (LRtol < N または LAtol < N の場合: スカラー, LRtol >= N かつ LAtol >= N の場合: 配列) Rtol は Atol と共に各ステップにおける局所誤差テストに用いられ, Y() の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. スカラーの場合 (i = 0 〜 N - 1): (Y(i) の局所誤差) <= Rtol(0)*Abs(Y(i)) + Atol(0) ただし, Rtol(0) と Atol(0) が同時に 0 であってはならない. 配列の場合 (i = 0 〜 N - 1): (Y(i) の局所誤差) <= Rtol(i)*Abs(Y(i)) + Atol(i) ただし, Rtol(i) と Atol(i) が同時に 0 であってはならない. [in] ATol() 配列 Atol(LAtol - 1) (LAtol >= 1) (Atol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値を表す. 本パラメータは LRtol および LAtol の値によりスカラーまたは配列を選択できる. (LRtol < N または LAtol < N の場合: スカラー, LRtol >= N かつ LAtol >= N の場合: 配列) Atol は Rtol と共に各ステップにおける局所誤差テストに用いられる (上記 Rtol 参照). [in] Mode 動作モード. 本プログラムは t0 から Tend までの積分を行うが, 中間結果の確認／出力の方法により4つの動作モードが提供される. どのモードでも T = Tend に達したときには積分を終了し, Info = 0 を返す. = 0: Tend まで戻らない. Tout は無視される. = 1: 成功したステップごとに戻る (Info = 1 を返す). Tout は無視される. T にはそのステップの終了点を返す. 最終ステップでは T = Tend となるようにステップ幅が調整される. = 2: 積分途中で Tout における中間結果を返すために戻る (T = Tout, Info = 2 を返す). 続けて次の Tout における解を求めるために積分を継続するためは, Tout を新たな値に変えて再度呼び出しを行うことができる. Tout 直前のステップでは T = Tout となるようにステップ幅が調節される. = 3: 積分途中で Tout における中間結果を返すために戻る (T = Tout, Info = 3 を返す). 続けて次の Tout における解を求めるために積分を継続するためは, Tout を新たな値に変えて再度呼び出しを行うことができる. Mode = 2 と異なり Tout における Y() の値は補間により求められ, Tout 直前のステップでのステップ幅の調節は行わず Mode = 1 と同じステップを踏む. [in,out] Info [in] 制御コード. = 0: 最初の呼び出し時(新たに問題を開始する場合)には Info = 0 と設定する. 全ての変数の初期化を行ってから計算を開始する. = 1, 2, 3: Info = 1, 2 または 3 で戻った場合, 計算を継続するために Tout だけを変更し, Info の値をそのままにして再呼び出しすることができる. [out] リターンコード. = 0: 正常終了. Tend までの積分が完了した. < 0: (-Info)番目の入力パラメータの誤り. = 1: Mode = 1 の中間結果出力のために戻った. 再呼び出しすることより次のステップまで進むことができる. = 2: Mode = 2 の中間結果出力のために T = Tout において戻った. Tout を再設定して再呼び出しすることができる. = 3: Mode = 3 の中間結果出力のために T = Tout において戻った. Tout を再設定して再呼び出しすることができる. = 11: (エラー) 最大ステップ数を超えた. = 12: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. = 13: (エラー) 行列が繰り返し特異になった. [out] Neval (省略可) 関数評価回数 (ヤコビ行列の計算のためのものは含まない). [out] Njac (省略可) ヤコビ行列評価回数 (有限差分の場合を含む). [out] Nstep (省略可) 全ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. [out] Ndec (省略可) 連立一次方程式のLU分解の回数. [out] Nsol (省略可) 連立一次方程式の解の計算回数. [in] Fjac (省略可) ヤコビ行列を求めるユーザー定義サブルーチンで, 次のように定義すること. (デフォルト値 = NullPtr) Sub Fjac(N As Long, T As Double, Y() As Double, Ypd() As Double) T および Y() におけるヤコビ行列を求め Ypd() に入れる. End Sub Ypd() (2次元配列) には, Mljac = N の場合にはフル行列形式(通常の n×n 2次元配列), Mljac < N の場合には帯行列形式(下記詳細参照)でヤコビ行列を格納すること. Fjac として NullPtr を指定すると, ユーザー定義サブルーチン呼び出しを行わず, ヤコビ行列を有限差分近似により求める. [in] Mljac (省略可) ヤコビ行列の下帯幅. (0 <= Mljac <= N) (デフォルト値 = N) Mljac = N の場合, ヤコビ行列はフル行列形式で格納される. Mljac < N の場合, ヤコビ行列は帯行列形式で格納される. [in] Mujac (省略可) ヤコビ行列の上帯幅. (0 <= Mujac <= N) (デフォルト値 = 0) Mljac = N の場合, Mujac は無視される. [in] Fmas (省略可) 質量マトリクス M を与えるユーザー定義サブルーチンで, 次のように定義すること. (デフォルト値 = NullPtr) Sub Fmas(N As Long, M() As Double) 質量マトリクスを M() に入れる. End Sub M() (2次元配列) には, Mlmas = N の場合にはフル行列形式, Mlmas < N の場合には帯行列形式で質量マトリクスを格納すること. Fmas として NullPtr を指定すると, ユーザー定義サブルーチン呼び出しを行わず, M は単位行列であるとみなす. [in] Mlmas (省略可) 質量マトリクス M の下帯幅. (0 <= Mljac <= N) (デフォルト値 = N) Mlmas = N の場合, M はフル行列形式で格納される. Mlmas < N の場合, M は帯行列形式で格納される. [in] Mumas (省略可) 質量マトリクス M の上帯幅. (0 <= Mumas <= N) (デフォルト値 = 0) Mlmas = N の場合, Mumas は無視される. [in] Hes (省略可) ヤコビ行列をヘッセンベルグ形に変換するかどうか指定する. 大規模でフル行列の場合にこの変換は有効である. 帯行列の場合および FMAS <> NullPtr の場合には無効である. = 0: 変換を行う. = 1: 変換を行わない. [in] Maxiter (省略可) 許される最大ステップ数. (デフォルト値 = 100000) [in] Nit1 (省略可) 陰的な方程式を解くために使用されるニュートン法の最大反復回数. (デフォルト値 = 7) [in] Startn (省略可) ニュートン法の出発値を指定する. (デフォルト値 = 0) = 0: 外挿した解を出発値として使う. = 1: 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 の変数がこの順番に現れるように作成されていなけらばならない. 誤差推定の際, 指数 2 変数には h を, 指数 3 変数には h^2 を乗算する. 常微分方程式 (ODE) の場合には Nind1 = N とする. [in] Pred (省略可) ステップ幅戦略を指定する. (デフォルト値 = 1) = 1: モデル予期コントローラ (グスタフソン). = 2: 伝統的なステップ幅コントロール. [in] M1,M2 (省略可) 連立微分方程式が次のような特殊な構造であるとする. y(i)' = y(i + m2) (i = 1, ..., m1) ただし, m1 は m2 の倍数であり, かつ残りの方程式が y'(m1), ..., y'(n - 1) に陽に依存しないものとする. このような場合, M1 = m1 (> 0) および M2 = m2 (> 0) (m1 + m2 <= n) と設定することにより計算時間を大幅に短縮することができる. (デフォルト値: M1 = 0, M2 = M1) [in] Nsmax,Nsmin (省略可) 最小段数 Nsmin および 最大段数 Nsmax (= 3, 5 または 7). (デフォルト値: Nsmin = 3, Nsmax = 7) 段数を ns とすると, その値により 2*ns - 1 次の陰的ルンゲ・クッタ法に対応する. [in] Ns (省略可) 最初のステップの ns の値. (3, 5 または 7) (デフォルト値 = Nsmin) [in] Cnt (省略可) カウンタ (Neval, Nstep, Naccept および Nreject) のリセット方法を指定する. (デフォルト値 = 0) = 0: Info = 0 のときにリセットする. = 1: Info = 0 であってもリセットしない. [in] Hinit (省略可) ステップ幅の初期値. (デフォルト値 = プログラムが自動推定する) [in] Hmax (省略可) ステップ幅の最大値. (デフォルト値 = Abs(Tend - T)) (Hmax = 0 であれば省略時の既定値とみなす) [in] Thet (省略可) ヤコビ行列を再計算するかどうかを指定する. (Thet < 1) (デフォルト値 = 0.001) ヤコビ行列の計算に時間がかかる場合には値を大きくすればよい (例えば 0.1). 小規模な場合には小さくすればよい (例えば 0.001). 負の値にするとステップが採用されるたびに強制的にヤコビ行列が再計算される. [in] Facl,Facr (省略可) ステップ幅選択パラメータ. (Facl <= 1, Facr >= 1) (デフォルト値: Facl = 0.2, Facr = 8) Facl <= hnew/hold <= Facr となるようにステップ幅が選ばれる. [in] Safe (省略可) ステップ幅推定時の安全係数. (0.001 < Safe < 1) (デフォルト値 = 0.9) [in] Quot1,Quot2 (省略可) Quot1 < hnew/hold < Quot2 であればステップ幅は変更されない. (Quot1 <= 1, Quot2 >= 1) (デフォルト値: Quot1 = 1, Quot2 = 1.2) 大規模な場合には, Quot1= 0.99, Quot2 = 2 が推奨される. これに加え Thet の値を大きくすることによりLU分解と計算時間を節約できる. 小規模な場合にはデフォルト値でよい. [in] Vitu (省略可) 縮小因子がこの値より小さいときに次数を増やす. (デフォルト値 = 0.002) [in] Vitd (省略可) 縮小因子がこの値より大きいときに次数を減らす. (デフォルト値 = 0.8) [in] Hhou,Hhod (省略可) ステップ幅の比が Hhou <= hnew/h <= Hhod を満たしたときにのみ次数を減らす. (デフォルト値: Hhou = 1.2, Hhod = 0.8) 文献 (1) 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)) 使用例 次の常微分方程式の初期値問題(スティフな問題)を解く. 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_Radaua() 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, Mode As Long, Neval As Long, Info As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Mode = 2 T = 0: Y(0) = 1: Y(1) = 2 Tend = 10 Info = 0 Do Tout = T + 1 Call Radaua(N, AddressOf F2, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Info, Neval) Debug.Print T, Y(0), Y(1) Loop While Info >= 1 And Info <= 3 Debug.Print Neval, Info End Sub 注 - 同じプログラムで Mode を 2 の代わりに 3 に変えれば密出力(補間)を使用する. 実行結果 1 0.367879441175671 0.908181747032157 2 0.135335283238606 -0.280811553265263 3 4.97870683691699E-02 -0.940205428205659 4 1.83156388907472E-02 -0.635327981830992 5 6.73794700079265E-03 0.290400132486203 6 2.47875217861059E-03 0.962649036304768 7 9.11881980026478E-04 0.75481411025709 8 3.3546262846512E-04 -0.145164571251474 9 1.23409804295512E-04 -0.91100685211013 10 4.53999298582475E-05 -0.839026129195407 270 0