Dehint Dehint_r Deoint Deoint_r Qagi Qagi_r Qawf Qawf_r Qk15i Qk15i_r

Dehint_r() Sub Dehint_r ( A As  Double, Result As  Double, Info As  Long, XX As  Double, YY As  Double, IRev As  Long, Optional Neval As  Long, Optional L As  Long, Optional Eps As  Double = -1  ) 半無限区間の積分 (自動積分) (二重指数関数(DE)公式) (リバースコミュニケーション版) 目的 本ルーチンは要求精度を満たす[a, +∞]におけるf(x)の積分値を求める. 被積分関数f(x)は変数IRevに従ってユーザーが必要な関数値を計算して与える. 積分値は二重指数関数(DE)公式を使用した自動積分により求められる. 二重指数関数(DE)公式(Definを参照)において以下の変換関数を使用することにより半無限区間[0, +∞]に対する公式を得る. (1) φ(t) = exp(t/2 - exp(-t)) (2) φ(t) = exp(t - exp(-t)) (3) φ(t) = exp(2sinh(t)) (1)は f(x) = f1(x)exp(-x^2) のように急激に減衰する場合に適している. (2)は f(x) = f2(x)exp(-x) のように指数的に減衰する場合に適している. (3)はゆるやかに減衰する有理関数や代数関数に適している. 本ルーチンでは被積分関数の挙動を調べて自動的に変換関数を選択して使用する(引数Lにどれを使用したか返す). 引数 [in] A 積分区間の下限 a. [out] Result 求められた積分値. [out] Info = 0: 正常終了. = -6: パラメータ IRev の誤り. = 1: 負側の減衰が遅いため収束判定条件を緩めた. = 2: 正側の減衰が遅いため収束判定条件を緩めた. = 3: 上の両方. = 4: 積分計算が収束しなかった. (Epsが小さすぎる, 区間内に特異点があるなどの原因による) [out] XX IRev = 1の場合, 関数値を求めるべき点を返す. [in] YY IRev = 1の場合, 再呼び出し時に関数値f(XX)を与えること. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックせよ. = 1: XXにおける関数値を求めYYに設定する. YY以外の変数を変更してはならない. [out] Neval (省略可) 被積分関数の評価回数. [out] L (省略可) 使用した変換関数. = 0: φ(t) = exp(t/2 - exp(-t)) = 1: φ(t) = exp(t - exp(-t)) = 2: φ(t) = exp(2sinh(t)) [in] Eps (省略可) 要求絶対誤差. (省略時 = 1.0e-14) (Eps <= 0 であれば省略時の既定値とみなす) 出典 森正武、「FORTRAN77数値計算プログラミング(増補版)」岩波書店 (1987) 使用例 次の定積分を求める. ∫ 1/(1 + x^2) dx [0, +∞] (= π/2) Sub Ex_Dehint_r() Dim A As Double, Result As Double, Info As Long Dim XX As Double, YY As Double, IRev As Long A = 0 IRev = 0 Do Call Dehint_r(A, Result, Info, XX, YY, IRev) If IRev = 1 Then YY = 1 / (1 + XX ^ 2) Loop While IRev <> 0 Debug.Print "S =", Result, "S(true) =", Dconst(13) / 2 Debug.Print "Info =", Info End Sub 実行結果 S = 1.5707963267949 S(true) = 1.5707963267949 Info = 0