◆ Defin()

Finite interval automatic quadrature (double exponential (DE) formula)

Purpose: This routine computes the integral of f(x) over [a, b], satisfying the requested accuracy, where f(x) is a given function defined by the user supplied routine F.
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]	F	The user supplied routine which computes the integrand function f(x) normally (when L = 0) defined as follows. Function F(X As Double) As Double F = computed f(X) value End Function When f(x) has singularities at endpoints, F may define the function f2(y) instead of f(x) and Defin may be called with L = 1 so that the integral of f(x) on [a, b] can be obtained with better precision, 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) X should not be changed in F.
[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]	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 defined routine F should define f2(y) instead of f(x). (See description of F) (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))

Function F1(X As Double) As Double

F1 = 1 / (1 + X ^ 2)

End Function

Sub Ex_Defin()

Dim A As Double, B As Double, Result As Double, Info As Long

A = 0: B = 4

Call Defin(AddressOf F1, A, B, Result, Info)

Debug.Print "S =", Result, "S(true) =", Atn(4)

Debug.Print "Info =", Info

End Sub

Defin
Sub Defin(F As LongPtr, A As Double, B As Double, Result As Double, Info As Long, Optional Neval As Long, Optional Eps As Double=-1, Optional L As Long=-1)
Finite interval automatic quadrature (double exponential (DE) formula)