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  ) Semi-infinite interval automatic quadrature for Fourier integrals (double exponential (DE) formula) (reverse communication version) Purpose This routine computes the following Fourier integrals by the automatic integration applying the double exponential (DE) formula. Ic = ∫ f(x) cosωx dx [a, +∞] or Is = ∫ f(x) sinωx dx [a, +∞] where f(x) is computed and provided by the user in accordance with IRev. Ic or Is is transformed to the intergral on the infinite interval [-∞, +∞] by using the following transformation function, and the integral value is computed by the trapezoidal rule. Ic: x = Mφ(t - π/2M)/ω Is: x = Mφ(t)/ω φ(t) can be selected from the following two functions by specifying the parameter. φ1(t) = t/(1 - exp(-2t -α(1 - exp(-t)) - β(exp(t) - 1))), β = 1/4, α = β/sqrt(1 + Mlog(1 + M)/4π) or φ2(t) = t/(1 - exp(-ksinh(t))), k = 2π The step size h and the constant M are selected to satisfy Mh = π. Please refer to the reference below for more details. Parameters [in] A Lower bound of integration interval a. [in] Omega ω value in cos or sin. [in] Integr Indicates which weight function is to be used. = 1: ∫ f(x) cosωx dx [0, +∞] = 2: ∫ f(x) sinωx dx [0, +∞] [out] Result Obtained approximation to the integral. [out] Info = 0: Successful exit = -3: The argument Integr had an illegal value (Integr <> 1 and Integr <> 2) = -8: The argument IRev had an illegal value. = 1: Not converged [out] XX When returned with IRev = 1, XX contains the abscissa where the function value should be evaluated and given in the next call. [in] YY When returned with IRev = 1, the function value f(XX) should be given in YY in the next call. [in,out] IRev Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered. [out] If IRev is not zero, complete the following tasks and call this routine again without changing IRev. = 0: Computation finished. See return code in Info. = 1: User should set the function values at XX in YY. Do not alter any variables other than YY. [out] Neval (Optional) Number of evaluations of f(x). [out] Err (Optional) Estimated relative error of obtained Result. To be overestimated for many cases. This routine may return as successful exit even if Err > Eps. [in] Eps (Optional) Relative accuracy requested. (default = 1.0e-10) (If Eps <= 0 or Eps > 1, the default value will be used) [in] Phi (Optional) The transformation function can be specified. (default = 0) = 0: φ1(t) = 1: φ2(t) (For other values, the default value will be used.) Reference 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. Example Program Compute the following integral. ∫ 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 Example Results S = 0.690194223521232 S(true) = 0.690194223521572 Info = 0