Debdf Debdf_r Radaua Radaua_r Rodasa Rodasa_r Seulexa Seulexa_r

Debdf_r() Sub Debdf_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, YYpd() As  Double, IRev As  Long, Optional Ijac As  Long = -1, Optional Ml As  Long = -1, Optional Mu As  Long = -1, Optional Mode As  Long = -1, Optional ITstop As  Long = -1, Optional Tstop As  Double  ) 常微分方程式の初期値問題 (1〜5次 後退微分公式 (BDF)) (リバースコミュニケーション版) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは微分方程式ソルバー・パッケージDEPACに収録されているプログラムのひとつである(DEPACは, Derkf, DeabmおよびDebdfからなる). Debdfは可変次数(1〜5次)後退微分公式のプログラムである. 3つの中で最も複雑である. 主に, スティフな微分方程式の低〜中精度な解を求めるために設計されている. 問題が強いスティフ性を持つ場合, DerkfおよびDeabmはDebdfに比べて極めて非効率である. しかし, 非スティフな問題に対しては, DebdfはDerkfおよびDeabmに比べて非効率である. これは, より多くのメモリを必要とし, よりオーバーヘッドが大きく, 次数が低いため効率的に精度を上げることができないためである. Debdfは A. C. Hindmarsh によるプログラム LSODE の修正版のドライバー・ルーチンである. Debdf_rはDebdfのリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. 新たなToutにおける解を求めるために積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, Tは前回のToutに等しくなければならない. また, Toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい. モードはパラメータModeで指定できる. [in] 独立変数Tの初期値. [out] 独立変数Tの最終ステップの最後の点の値 (インターバル・モードでは通常Toutに等しい). 解がこの点まで正常に求められたことを示す. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数Y()の初期値. [out] 最終のT(インターバル・モードにおいてはToutに等しい)において求められた解(の近似値). [in] Tout 解を求めたい点を設定する. Tout = T とすることができ, その場合はTにおける微分係数を計算して戻る. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても本来のとおりToutを指定しなければならない. [in] RTol RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)または配列(LRTol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. [in] ATol ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)または配列(LATol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. [in,out] Info [in] 制御コード. = 0: 問題の最初の呼び出し時(新たに問題を開始する場合)には Info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される. = 1, 2, 11〜15: Info = 1, 2, 11〜15 で戻った場合, 計算を継続するために, Infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ). [out] リターンコード. ユーザーはこの値をチェックし, Info = 1, 2, 11〜15 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる. = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ T の誤り. (継続呼び出しで T = Tout または T <> 前回のTout) = -4: パラメータ Y() の誤り. = -5: パラメータ Tout の誤り. (積分の方向を変更した) = -6: パラメータ Rtol() の誤り. (RTol(i) < 0) = -7: パラメータ Atol() の誤り. (ATol(i) < 0) = -8: パラメータ Info の誤り. (Info <> 0, 1, 2, 11〜15) = -10: パラメータ MlまたはMu の誤り. (Ml < 0, Ml >= N, Mu < 0 または Mu >= N) = 1: 正常終了 (T = Tout). Toutを再設定して再呼び出しすることができる. ToutはTと異なる値であること. = 2: 中間結果出力モードのため戻った (まだ Toutに達していない). Tout方向の次のステップを続行するために再呼び出しすることができる. = 11: 最大ステップ数(10000)を超えた. 再呼び出しして続行することができる. さらに10000ステップ続けることができる. = 12: 誤差許容値が厳しすぎる. 許容値を変更してから再呼び出しすることができる. = 13: 求められた解が0であるため, 純粋な相対誤差テスト(ATol = 0)ができない. Atolを正の値にして再呼び出しして続行することができる. = 14: 最後の試行ステップにおいて収束テストの失敗を繰り返した. Infoを変えずに再呼び出しして続行することができるが, ヤコビ行列が正しくないことが原因かもしれない. = 15: 最後の試行ステップにおいて誤差テストの失敗を繰り返した. Infoを変えずに再呼び出しして続行することができるが, 解に特異性があることが原因かもしれない. = 16: 無限ループが検出された. [out] TT IRev = 1, 2: 再呼び出し時に微分係数値あるいはヤコビ行列を求めてYYp()あるいはYYpd()に入れるべき点のTの値を返す. [out] YY() 配列 YY(LYY - 1) (LYY >= N) IRev = 1, 2: 再呼び出し時に微分係数値あるいはヤコビ行列を求めてYYp()あるいはYYpd()に入れるべき点の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 >= N (フル行列のとき) または Ml+Mu+1 (帯行列のとき), LYYpd2 >= N) IRev = 2: T(= TT)およびY(= YY())におけるヤコビ行列(dfi/dyj)を計算して, 再呼び出し時にに入力する. Ml = N の場合 YYpd()は通常のフル行列, 0 <= Ml < N の場合帯行列形式である. [in,out] IRev リバースコミュニケーションの制御変数 [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: TTおよびYY()における微分係数の計算値をYYp()に設定する. YYp()以外の変数を変更してはならない. = 2: TTおよびYY()におけるヤコビ行列の計算値をYYpd()に設定する. YYpd()以外の変数を変更してはならない. [in] Ijac (省略可) ヤコビ行列の計算方法. (省略時 = 0) = 0: 差分近似する. (IRev = 2では戻らない) = 1: 計算が必要なときにIRev = 2で戻り, ユーザーが計算して入力する. (他の値であれば省略時の既定値とみなす.) [in] Ml (省略可) Ml = N の場合, ヤコビ行列はN×Nフル(通常の)行列として扱われる. 0 <= Ml < N であればヤコビ行列を帯行列形式で扱う(2*Ml + Mu < N のときに効果がある). その場合, Mlは下帯幅を表す. (Ml >= 0, Ml <= N) (省略時 = N) [in] Mu (省略可) ヤコビ行列を帯行列として扱う場合に上帯幅を表す. Ml = N の場合には参照されない. (Mu >= 0, Mu < N) (省略時 = 0) [in] Mode (省略可) 動作モード. (省略時 = 0) = 0: toutで戻る (インターバル・モード). = 1: ステップごとに戻る (中間結果出力モード) (他の値であれば省略時の既定値とみなす.) [in] ITstop (省略可) 本ルーチンはtoutを超えて積分を行い内挿によりtoutにおける値を求めることがあるが, ある点において不連続になったり, その点を超えて関数あるいは微分が定義されていなかったりして, その点を超えて積分を行うのが許されないことがある. そのような場合, ストップポイント Tstop を設定することができ, その点を超えての積分を行わないようにすることができる. (省略時 = 0) = 0: ストップポイントを設定しない = 1: Tstopをストップポイントとして設定する (上記以外の値であれば省略時の既定値とみなす) [in] Tstop (ITstop = 0 のとき省略可) ストップポイントの値. 出典 SLATEC (DEPAC) 使用例 次の常微分方程式の初期値問題(スティフな問題)を解く. 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_Debdf_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, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, YYpd(0) As Double, 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 Debdf_r(N, T, Y(), Tout, RTol(), ATol(), Info, TT, YY(), YYp(), YYpd(), IRev) If IRev = 1 Then Call F2(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Derkf_r: Info =", Info Exit Sub End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879442394948 0.908181748259189 2 0.135335284386392 -0.280811552155648 3 0.049787069049436 -0.940205427530908 4 1.83156392354596E-02 -0.635327981611744 5 6.73794718964386E-03 0.290400132587821 6 2.47875225773856E-03 0.96264903888483 7 9.11881999298641E-04 0.75481413632096 8 3.35462642293389E-04 -0.145164571183809 9 1.2340981022923E-04 -0.911006852067105 10 4.53999332058006E-05 -0.839026129132657