defint defint_r qag qag_r qags qags_r qng qng_r

defint_r() void defint_r ( double  a, double  b, double  eps, int  l, double *  result, int *  neval, int *  info, double *  xx, double  yy, int *  irev  ) 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. [in] b Upper limit of integration. [in] eps Absolute accuracy requested. max(|eps|, 1.0e-32) is used as the tolerance. [in] l = 0: yy returnes the computed value of f(xx). = 1: yy returnes the computed value of f2(xx). (For other values, l = 0 is assumed) [out] result Approximation to integral of f(x) over [a, b]. [out] neval Number of integrand evaluations. [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 [out] xx irev = 1 to 5: xx contains the abscissa where the function value should be evaluated and given in the next call. [in] yy irev = 1 to 5: The function value f(xx) (when l = 0) or f2(xx) (when l = 1) should be given in yy in the next call, where 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) The latter option (l = 1) may be used when f(x) has singularities at endpoints. User may provide the function value of f2(y) instead of f(x) and call defint with l = 1 so that the integral of f(x) on [a, b] can be obtained with better precision. [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. = 1 to 5: User should set the function value at xx in yy. Do not alter any variables other than yy. Reference Masatake Mori "FORTRAN77 Numerical Calculation Programming (augmented edition)" Iwanami Shoten, 1987. (Japanese book)