Deabm Deabm_r Derkfa Derkfa_r Dop853a Dop853a_r Dopn1210 Dopn1210_r Dopn43 Dopn43_r Dopn64 Dopn64_r Dopn86 Dopn86_r Dopri5a Dopri5a_r Dverka Dverka_r Odex2a Odex2a_r Odexa Odexa_r Retarda Retarda_r Ylaga

Deabm() Sub Deabm ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Mode As  Long = -1, Optional ITstop As  Long = -1, Optional Tstop As  Double  ) 常微分方程式の初期値問題 (1〜12可変次数 アダムス・バシュフォース・ムルトン法) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは微分方程式ソルバー・パッケージDEPACに収録されているプログラムのひとつである(DEPACは, Derkf, DeabmおよびDebdfからなる). Deabmは可変次数(1〜12次)アダムス・バシュフォース・モルトンの予測子・修正子法のプログラムである. 複雑さはDerkfとDebdfの中間ぐらいである. 主に, 微分係数の計算に時間がかかり, 高精度な解を求めたい, あるいは, 非常に多くの点における解を求めたい非スティフおよびややスティフな微分方程式を解くために設計されている. Deabmは L. F. Shampine および M. K. Gordon によるプログラム ODE の修正版のドライバー・ルーチンである. 引数 [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()以外の変数を変更してはならない. [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〜14: Info = 1, 2, 11〜14 で戻った場合, 計算を継続するために, Infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ). [out] リターンコード. ユーザーはこの値をチェックし, Info = 1, 2, 11〜14 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる. = -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〜14) = 1: 正常終了 (T = Tout). Toutを再設定して再呼び出しすることができる. ToutはTと異なる値であること. = 2: 中間結果出力モードのため戻った (まだ Toutに達していない). Tout方向の次のステップを続行するために再呼び出しすることができる. = 11: 最大ステップ数(10000)を超えた. 再呼び出しして続行することができる. さらに10000ステップ続けることができる. = 12: 誤差許容値が厳しすぎる. 許容値を変更してから再呼び出しすることができる. = 13: 求められた解が0であるため, 純粋な相対誤差テスト(ATol = 0)ができない. Atolを正の値にして再呼び出しして続行することができる. = 14: スティフな問題の可能性がある. 再呼び出しして続行することができる. しかし, スティフな問題を効率的に扱うことができるDEPACのDebdfのような他のプログラムが推奨される. = 16: 無限ループが検出された. [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 = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T) End Sub Sub Ex_Deabm() 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 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 Call Deabm(N, AddressOf F1, T, Y(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Deabm: Info =", Info Exit Sub End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441163166 0.908181747036676 2 0.135335283217889 -0.280811553319014 3 4.97870683619061E-02 -0.940205428247124 4 1.83156388898816E-02 -0.63532798198588 5 6.73794699179997E-03 0.29040013250825 6 2.47875218622002E-03 0.962649038865161 7 9.11881990311618E-04 0.754814136279958 8 3.35462640476657E-04 -0.145164571214819 9 1.23409802287714E-04 -0.911006852099835 10 4.53999261671108E-05 -0.839026129149036