Defin Defin_r Qag Qag_r Qags Qags_r Qng Qng_r

Defin_r() Sub Defin_r ( A As  Double, B As  Double, Result As  Double, Info As  Long, XX As  Double, YY As  Double, IRev As  Long, Optional Neval As  Long, Optional Eps As  Double = -1, Optional L As  Long = -1  ) Finite interval automatic quadrature (double exponential (DE) formula) (reverse communication version) Purpose This routine computes the integral of f(x) over [a, b], satisfying the requested accuracy, where f(x) is a given function. User should provide the necessary computed values of f(x) according to the argument IRev. The result is obtained by the automatic integration applying the double exponential (DE) formula. The integral interval [-1, 1] can be transformed into the infinite interval [-∞, ∞] by the variable transformation. ∫ f(x)dx [-1, 1] = ∫ g(t)dt [-∞, ∞] where x=φ(t), g(t) = f(φ(t))φ'(t) φ(t) is the monotone function, and it satisfies φ(-∞) = -1, φ(∞) = 1. By choosing suitable φ(t), g(t) can have the double exponential asymptotic behavior and can be integrated efficiently by trapezoidal rule. The following transformation function is used in this routine. φ(t) = tanh((π/2)sinh(t)) The DE formula can be applied even if the singularities exist at the end points. Parameters [in] A Lower limit of integration a. [in] B Upper limit of integration b. [out] Result Approximation to the integral. [out] Info = 0: Successful exit. = 1: Slow decay on negative side. = 2: Slow decay on positive side. = 3: Both of above. = 4: Insufficient mesh refinement. (Eps too small, singularity within integral interval, etc.) [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 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. Normally (L = 0) set f(XX) as the function value in YY. If f(x) has singularities at endpoints, set L = 1 and set f2(XX) as the function value in YY so that the integral of f(x) can be obtained with better precision. f2(y) is defined as follows.   f2(y) = f(a - y) (-(b - a)/2 < y < 0)   f2(y) = f(b - y) (0 < y <= (b - a)/2) [out] Neval (Optional) Number of integrand evaluations. [in] Eps (Optional) Absolute accuracy requested. (default = 1.0e-14) (If Eps <= 0, the default value will be used) [in] L (Optional) Definition of F. (default = 0) = 0: Normal. = 1: The integrand function will be integrated after transformation of variable to avoid the computation on the limit. The user should compute f2(y) instead of f(x). (See description of IRev) (For other values, the default value is assumed) Reference (Japanese book) Masatake Mori "FORTRAN77 Numerical Calculation Programming (augmented edition)" Iwanami Shoten (1987) Example Program Compute the following integral. ∫ 1/(1 + x^2) dx [0, 4] (= atan(4)) Sub Ex_Defin_r() Dim A As Double, B As Double, Result As Double, Info As Long Dim XX As Double, YY As Double, IRev As Long A = 0: B = 4 IRev = 0 Do Call Qag_r(A, B, Result, Info, XX, YY, IRev) If IRev = 1 Then YY = 1 / (1 + XX ^ 2) Loop While IRev <> 0 Debug.Print "S =", Result, "S(true) =", Atn(4) Debug.Print "Info =", Info End Sub Example Results S = 1.32581766366803 S(true) = 1.32581766366803 Info = 0