Dehint Dehint_r Deoint Deoint_r Qagi Qagi_r Qawf Qawf_r Qk15i Qk15i_r

Deoint() Sub Deoint ( F As  LongPtr, A As  Double, Omega As  Double, Integr As  Long, Result As  Double, Info 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)はユーザー定義サブルーチンFにより与えられる関数である. 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] F 被積分関数f(x)を求めるユーザー定義サブルーチンで, 次のように定義すること. Function F(X As Double) As Double F = f(X) End Function Xを変更しないこと. [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: 正常終了 = -4: パラメータ Integr の誤り (Integr <> 1 かつ Integr <> 2) = 1: 収束しなかった [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*√π) Function F1(X As Double) As Double F1 = Exp(-X ^ 2) End Function Sub Ex_Deoint() Dim A As Double, Omega As Double, Integr As Long, Result As Double Dim Info As Long, Neval As Long A = 0 Omega = 1 Integr = 1 Call Deoint(AddressOf F1, A, Omega, Integr, Result, Info) 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