Dehint Dehint_r Deoint Deoint_r Qagi Qagi_r Qawf Qawf_r Qk15i Qk15i_r

Deoint_r() Sub Deoint_r ( A As  Double, Omega As  Double, Integr As  Long, Result As  Double, Info As  Long, XX As  Double, YY As  Double, IRev As  Long, Optional Neval As  Long, Optional Err As  Double, Optional Eps As  Double = -1, Optional Phi As  Long = 0  ) 半無限区間の積分 (フーリエ型積分) (自動積分) (二重指数関数(DE)公式) (リバースコミュニケーション版) 目的 本ルーチンは二重指数関数(DE)公式を使用した自動積分によりフーリエ型積分を求める. Ic = ∫ f(x) cosωx dx [a, +∞] または Is = ∫ f(x) sinωx dx [a, +∞] ここで, f(x)は変数IRevに従ってユーザーが必要な関数値を計算して与える. DE公式では次の変換関数を使用することにより無限区間[-∞, +∞]の積分に変換して台形公式により積分値を計算する. Ic: x = Mφ(t - π/2M)/ω Is: x = Mφ(t)/ω φ(t)はパラメータの指定により次の2種類から選択することができる. φ1(t) = t/(1 - exp(-2t -α(1 - exp(-t)) - β(exp(t) - 1))), β = 1/4, α = β/sqrt(1 + Mlog(1 + M)/4π) または φ2(t) = t/(1 - exp(-ksinh(t))), k = 2π 刻み幅hと定数Mは, Mh = π が成り立つように選ばれる. 詳細は下記文献を参照のこと. 引数 [in] A 積分の下限値 a. [in] Omega cosまたはsinの中のωの値. [in] Integr 重み関数を指定する. = 1: ∫ f(x) cosωx dx [0, +∞] = 2: ∫ f(x) sinωx dx [0, +∞] [out] Result 求められた積分値. [out] Info = 0: 正常終了 = -3: パラメータ Integr の誤り (Integr <> 1 かつ Integr <> 2) = -8: パラメータ IRev の誤り. = 1: 収束しなかった [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 (省略可) f(x)の評価回数. [out] Err (省略可) Resultの推定相対誤差. 多くの場合大きめに算出される. また, Err > Eps であっても正常終了として戻ることがある. [in] Eps (省略可) 要求相対誤差. (省略時 = 1.0e-10) (Eps <= 0 または Eps > 1 であれば省略時の既定値とみなす) [in] Phi (省略可) 変換関数を指定する. (省略時 = 0) = 0: φ1(t) = 1: φ2(t) (上記以外の値であれば省略時の既定値とみなす) 参考文献 Ooura and Mori, "The double exponential formula for oscillatory functions over the half infinite interval", J. Comput. Appl. Math., 38 (1991), 353-360. Ooura and Mori, "A robust double exponential formula for Fourier-type integrals", J. Comput. Appl. Math., 112 (1999), 229-241. 使用例 次の定積分を求める. ∫ exp(-x^2)cos(x) dx [0, +∞] (= exp(-1/4)*1/2*√π) Sub Ex_Deoint_r() Dim A As Double, Omega As Double, Integr As Long, Result As Double, Info As Long Dim XX As Double, YY As Double, IRev As Long A = 0 Omega = 1 Integr = 1 IRev = 0 Do Call Deoint_r(A, Omega, Integr, Result, Info, XX, YY, IRev) If IRev = 1 Then YY = Exp(-XX ^ 2) Loop While IRev <> 0 Debug.Print "S =", Result, "S(true) =", Exp(-1 / 4) * Dconst(17) / 2 Debug.Print "Info =", Info End Sub 実行結果 S = 0.690194223521232 S(true) = 0.690194223521572 Info = 0