NuMFor 9f2ab49 (2024-04-08)
Numerical (Modern) Fortran. Library for Simple Numerical computing
Submodule Integrate

There are implemented several routines for integration of functions, and also sampled values via Trapezoid and Simpson methods.

Each algorithm computes an approximation to a definite integral of the form,

\[J \approx I := \int_a^b f(x) w(x) dx\]

where $w(x)$ is a weight function, which may be $w(x) = 1$.

Integrable functions

All routines described here may integrate functions that return either real or complex values. They may have one of the following signatures:

  • real function f(x) of a real argument x
  • complex function f(x) of a real argument x
  • real function f(x,args) of real arguments x, and array args
  • complex function f(x,args) of real arguments x, and array args

For functions that need extra arguments besides the integration variable x, the argument array must be passed to the integration routines.

There is also the possibility of using

Trapezoid and Simpson rules

There are implemented routines for integration of functions by simple Trapezoid and Simpson rules via two routines trapz and simps. The use of both routines is very similar. For instance:

If we have data sampled to a grid of ordered values, we can integrate them by using trapezoid or simpson quadrature similarly as with functions. In order to use them we have only to call trapz or simps

program ex_trapz
USE numfor, only: dp, zero, m_pi, str, trapz, linspace
implicit none
real(dp), dimension(20) :: x
real(dp), dimension(20) :: y
real(dp) :: Integ1
intrinsic dsin
integ1 = trapz(dsin, zero, m_pi)
print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
! \int \sin(x) dx = 1.999832163894 (Difference=0.000167836106)
! Integration of sampled data
x = linspace(zero, m_pi, 20)
y = dsin(x)
integ1 = trapz(y, zero, m_pi)
print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
! \int \sin(x) dx = 1.9980436909706 (Difference=0.0019563090294)
! Integration of a function
integ1 = trapz(dsin, zero, m_pi, n=size(y))
print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
! \int \sin(x) dx = 1.9980436909706 (Difference=0.0019563090294)
end program ex_trapz

The use of Simpson routine is similar, and need only to replace trapz() by simps()

Adaptive Simpson method

In addition to The simple, fixed points, Simpson quadrature routines there are implemented Globally Adaptive Simpson routines for finite and infinite integration domains.

Routine iads follows a scheme in which the integration interval is bisected, and the Simpson integration result for the sum of the subintervals is compared with the result for the parent interval. For each subdivision, only three additional function evaluations are needed. This routine is regarded as a good choice when moderate precision ( $10^{-3}, 10^{-4}$) is required.

program ex_iads
USE numfor, only: dp, zero, m_pi, str, iads
implicit none
real(dp) :: err
integer :: N
real(dp) :: Integ1
intrinsic dsin
call iads(dsin, zero, m_pi, integ1)
print "(A)", 'integrate(sin(x), 0, pi) = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
! integrate(sin(x), 0, pi) = 1.9999998386284 (Difference=1.6137164227104e-07)
! < [Simple]
call iads(fquad451, 1.e-12_dp, m_pi, integ1, epsrel=1.e-7_dp, abserr=err, neval=n)
print "(A)", 'integrate( \log(x) \sqrt(x), 1.e-12, pi) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(abs(n))
! integrate( \log(x) \sqrt(x), 1.e-12, pi) = 1.7746775593141 (Error=1.3697319985965e-05) with N = 89
contains
function fquad451(x) result(y)
implicit none
real(dp) :: y
real(dp), intent(IN) :: x
y = sqrt(x) * log(x)
end function fquad451
end program ex_iads

Routine iadsi is available for semi-infinite integration domains.

For instance, for the function

we see that the main contribution comes from the region below $x \sim 1$. Then, we can use that information as argument brkpts, as shown in the following example

program ex_iadsi
USE numfor, only: dp, zero, m_pi, str, iadsi, nf_inf, qags
implicit none
real(dp) :: err
integer :: N
real(dp) :: Integ1
intrinsic dsin
call iadsi(f, zero, [1._dp], integ1, epsrel=1.e-7_dp, abserr=err, neval=n)
print "(A)", 'integrate(log(1+x)/(1 + 100*x^2), 0, inf) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(abs(n))
! integrate(log(1+x)/(1 + 100*x^2), 0, inf) = 0.0335216255386 (Error=1.259162853146e-12) with N = 165
contains
function f(x) result(y)
implicit none
real(dp) :: y
real(dp), intent(IN) :: x
y = log(1 + x) / (1 + 100 * x**2)
end function f
end program ex_iadsi

Integration using the Tanh-sinh scheme

Routine qnthsh is an implementation of the Tanh-sinh quadrature method. As explained in Wikipedia, the tanh-sinh or Double exponential quadrature is a method for numerical integration introduced by Hidetosi Takahasi and Masatake Mori in 1974. It is based in the change of variables

\[ x := g(t) = \tanh\left(\frac{\pi}{2}\sinh t\right)\,\]

to transform an integral on the interval x ∈ (−1, +1) to an integral on the entire real line t ∈ (−∞, +∞). After this transformation, the integrand decays with a double exponential rate, and thus, this method is also known as the "Double Exponential (DE) formula" (Mori, Masatake (2005), "Discovery of the Double Exponential Transformation and Its Developments", Publications of the Research Institute for Mathematical Sciences, 41 (4): 897-935).

For a given step size h, the integral is approximated by the sum

\[ \int_{-1}^{1} f(x) \, dx \approx \sum_{k=-\infty}^\infty w_{k} f(x_{k}),\]

with the abscissas

\[ x_k := g(h k) = \tanh\left(\frac{\pi}{2}\sinh kh\right) \]

and the weights

\[w_k := g'( h k) = \frac{\frac{\pi}{2}h \cosh(kh)}{\cosh^{2} \big(\frac{\pi}{2} \sinh (kh) \big)}.\]

The Tanh-Sinh method is quite insensitive to endpoint behavior. Should singularities or infinite derivatives exist at one or both endpoints of the (−1, +1) interval, these are mapped to the (−∞,+∞) endpoints of the transformed interval, forcing the endpoint singularities and infinite derivatives to vanish. This results in a great enhancement of the accuracy of the numerical integration procedure, which is typically performed by the Trapezoidal Rule. In most cases, the transformed integrand displays a rapid roll-off (decay), enabling the numerical integrator to quickly achieve convergence.

Examples of use

The most simple example of using the routine is given by

1program ex_qthsh1
2 USE numfor, only: dp, zero, m_pi, str, center, qnthsh
3
4 real(dp) :: err
5 integer :: N
6
7 real(dp) :: Integ1
8 intrinsic dsin
9 print "(A)", center(" Integrate sin(x) between 0 and pi ", 70, '-')
10 call qnthsh(dsin, zero, m_pi, integ1)
11 print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
12 ! that prints:
13 !------------------ Integrate sin(x) between 0 and pi -----------------
14 ! \int \sin(x) dx = 2 (Difference=0)
15
16 print "(A)", center(" Get estimation of error and number of evaluations ", 70, '-')
17 call qnthsh(dsin, zero, m_pi, integ1, abserr=err, neval=n)
18 print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
19 ! that prints:
20 !---------- Get estimation of error and number of evaluations ---------
21 ! \int \sin(x) dx = 2 (Error=4.8361314952672e-12) with N = 65
22end program ex_qthsh1

Here in the call in line 10 is the simplest possible and only uses the required arguments.

Quadpack routines

The code implemented here has been slightly adapted from the original Quadpack routines (by Piessens, de Doncker-Kapenga, Ueberhuber and Kahaner). It has been ported to Modern Fortran by following several sources: Gnu Scientific Library, J. Burkardt's page, and SciFortran. Moreover, the routines have been tweaked to work with both real and complex functions.

There are routines for adaptive and non-adaptive integration of general functions, with specific routines for some particular difficult cases. These include integration over infinite and semi-infinite ranges, singular integrals, including potential and logarithmic singularities, computation of Cauchy principal values and oscillatory integrals.

The algorithms in Quadpack use a naming convention based on the following letters:

Symbol Explanation
Q quadrature routine
N non-adaptive integrator
A adaptive integrator
G general integrand (user-defined)
W weight function with integrand
S singularities can be more readily integrated
P points of special difficulty can be supplied
I infinite range of integration
O oscillatory weight function, cos or sin
F Fourier integral
C Cauchy principal value

The algorithms are built on pairs of quadrature rules, a higher order rule and a lower order rule. The higher order rule is used to compute the best approximation to an integral over a small range. The difference between the results of the higher order rule and the lower order rule gives an estimate of the error in the approximation.

All subroutines present a similar interface. They accept several optional input and output arguments but, in their simpler form, all integrators without weight functions may be called as

call integrator(f, a, b, IntVal)

where f is a real or complex function f(x,[args]) as described in Integrable functions, IntVal is the value returned with an approximation to the integral of f over the (possibly infinite) interval (a,b).

There are available the following integration routines:

Subroutine QNG for non-adaptive integration

The subroutine qng is a simple non-adaptive automatic integrator, based on a sequence of rules with increasing degree of algebraic precision (Patterson, 1968). It applies the Gauss-Kronrod 10-point, 21-point, 43-point and 87-point integration rules in succession until an estimate of the integral is achieved within the desired absolute and relative error limits. Each successive approximation reuses the previous function evaluations.

Examples of use

1program ex_qng
2 USE numfor, only: dp, zero, m_pi, str, center, qng
3 implicit none
4
5 real(dp) :: err
6 integer :: N
7 real(dp) :: w
8
9 real(dp) :: Integ1
10 complex(dp) :: Integ2
11 intrinsic dsin
12 print "(A)", center(" Integrate sin(x) between 0 and pi ", 70, '-')
13 call qng(dsin, zero, m_pi, integ1)
14 print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
15 ! that prints:
16 !------------------ Integrate sin(x) between 0 and pi -----------------
17 ! \int \sin(x) dx = 2 (Difference=0)
18
19 w = 100._dp
20 print "(A)", center(" Get estimation of error and number of evaluations ", 70, '-')
21 call qng(fquad451, zero, m_pi, integ1, abserr=err, neval=n)
22 print "(A)", '\int \sqrt(x) \log(x) dx = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
23 ! that prints:
24 !---------- Get estimation of error and number of evaluations ---------
25 ! \int \sqrt(x) \log(x) dx = 1.774674495015 (Error=6.9889313343983e-05) with N = 87
26 print "(A)", center(" Integrate a complex function ", 70, '-')
27 call qng(fquad452wc, zero, m_pi, [w], integ2, abserr=err, neval=n)
28 print "(A)", '\int \cos(w \sin(x)) dx = '//str(integ2)//" (Error="//str(err)//") with N = "//str(n)
29 ! that prints:
30 !-------------------- Integrate a complex function --------------------
31 ! \int \cos(w \sin(x)) dx = (0.062787400402-0.2226721656122j) (Error=2.8652541230259) with N = 87
32
33contains
35 function fquad451(x) result(y)
36 implicit none
37 real(dp) :: y
38 real(dp), intent(IN) :: x
39 y = sqrt(x) * log(x)
40 end function fquad451
41
43 function fquad452wc(x, w) result(y)
44 implicit none
45 complex(dp) :: y
46 real(dp), intent(IN) :: x
47 real(dp), dimension(:), intent(IN) :: w
48 y = cmplx(cos(w(1) * sin(x)), sin(w(1) * sin(x)), kind=dp)
49 end function fquad452wc
50
51end program ex_qng

Subroutine QAG for simple adaptive integration

The subroutine qag uses a global adaptive scheme in which the interval is divided and at each step the integral is calculated for each subinterval, and the one with the larger error is bisected. This scheme is recusively applied and the most difficult points are calculated in finer grids.

Examples of use

1program ex_qag
2 USE numfor, only: dp, zero, m_pi, str, center, qag
3 implicit none
4
5 real(dp) :: err
6 integer :: N
7 real(dp) :: w
8
9 real(dp) :: Integ1
10 complex(dp) :: Integ2
11 real(dp), parameter :: epsrel = 1.e-3_8
12
13 intrinsic dsin
14 print "(A)", center(" Integrate sin(x) between 0 and pi ", 70, '-')
15 call qag(dsin, zero, m_pi, integ1)
16 print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
17 ! that prints:
18 !------------------ Integrate sin(x) between 0 and pi -----------------
19 ! \int \sin(x) dx = 2 (Difference=0)
20
21 w = 100._dp
22 print "(A)", center(" Get estimation of error and number of evaluations ", 70, '-')
23 call qag(fquad451, zero, m_pi, integ1, abserr=err, neval=n)
24 print "(A)", '\int \sqrt(x) \log(x) dx = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
25
26 call qag(fquad451, zero, m_pi, integ1, rule='qk21', abserr=err, neval=n)
27 print "(A)", '\int \sqrt(x) \log(x) dx = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
28 ! ---------- Get estimation of error and number of evaluations ---------
29 ! \int \sqrt(x) \log(x) dx = 1.7746752010439 (Error=7.9254688067762e-06) with N = 375
30 ! \int \sqrt(x) \log(x) dx = 1.7746751858103 (Error=1.4425452046185e-05) with N = 441
31
32 print "(A)", center(" Integrate a complex function ", 70, '-')
33 call qag(fquad452wc, zero, m_pi, [w], integ2, rule='qk61', abserr=err, neval=n)
34 print "(A)", '\int \cos(w \sin(x)) dx = '//str(integ2)//" (Error="//str(err)//") with N = "//str(n)
35 ! that prints:
36 !-------------------- Integrate a complex function --------------------
37 ! \int \cos(w \sin(x)) dx = (0.062787390876-0.2226721893626j) (Error=1.7375840100087e-08) with N = 427
38
39contains
41 function fquad451(x) result(y)
42 implicit none
43 real(dp) :: y
44 real(dp), intent(IN) :: x
45 y = sqrt(x) * log(x)
46 end function fquad451
47
49 function fquad452wc(x, w) result(y)
50 implicit none
51 complex(dp) :: y
52 real(dp), intent(IN) :: x
53 real(dp), dimension(:), intent(IN) :: w
54 y = cmplx(cos(w(1) * sin(x)), sin(w(1) * sin(x)), kind=dp)
55 end function fquad452wc
56
57end program ex_qag

Subroutine QAGS for general integration with singularities

Subroutine qags is a global adaptive algorithm with bisection of the integration interval and a Wynn-epsilon algorithm to extrapolate and speed-up the convergence of the integral. This routine is able to integrate functions with singularities of many types.

The limits of integration may be finite or infinite. By default it uses a Gauss-Kronrod rule of 21 points for finite limits, and a 15-points rule for infinite limits. The rule may be given to the routine as an argument.

The integration over the infinite intervals are is performed by means of mappings to the semi-open interval $(0,1]$.

Examples of use

1program ex_qags
2 USE numfor, only: dp, zero, m_pi, nf_inf, nf_minf, str, center, qags
3 implicit none
4
5 real(dp) :: err
6 integer :: N
7 real(dp) :: w
8
9 real(dp) :: Integ1
10
11 intrinsic dsin
12 print "(A)", center(" Integrate sin(x) between 0 and pi ", 70, '-')
13 call qags(dsin, zero, m_pi, integ1)
14 print "(A)", 'integrate(sin(x), 0, pi) = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
15 ! that prints:
16 !------------------ Integrate sin(x) between 0 and pi -----------------
17 ! integrate(sin(x), 0, pi) = 1.9999999827467 (Difference=1.725330678326e-08)
18
19 w = 100._dp
20
21 ! Get estimation of error and number of evaluations
22 call qags(fquad453, zero, 1._dp, integ1, gkrule='qk15', abserr=err, neval=n)
23 print "(A)", 'integrate( \log(x)/\sqrt(x), 0, 1) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
24 ! integrate( \log(x)/\sqrt(x), 0, 1) = -3.9999999999999 (Error=5.1247894816697e-13) with N = 225
25
26 call qags(fquad455w, zero, nf_inf, [w], integ1, gkrule='qk21', abserr=err, neval=n)
27 print "(A)", ∞'integrate( \log(x)/(1 + w x^2), 0, +) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
28 ! integrate( \log(x)/(1 + w x^2), 0, +∞) = -0.3616892188416 (Error=2.170181735946e-07) with N = 399
29contains
30
32 function fquad453(x) result(y)
33 implicit none
34 real(dp) :: y
35 real(dp), intent(IN) :: x
36 y = log(x) / sqrt(x)
37 end function fquad453
38
40 function fquad455w(x, w) result(y)
41 implicit none
42 real(dp) :: y
43 real(dp), intent(IN) :: x
44 real(dp), dimension(:), intent(IN) :: w
45 y = log(x) / (1 + w(1) * x**2)
46 end function fquad455w
47end program ex_qags

Subroutine QAGP integration with known singular points

Subroutine qagp uses the same strategy than qags to integrate a function that presents singularities in the interior of the integration domain. The user has to supply an array with the points where difficulties arise.

Examples of use

1program ex_qagp
2 USE numfor, only: dp, zero, str, qagp
3 implicit none
4
5 real(dp) :: err
6 integer :: N
7 real(dp), parameter :: epsabs = 0.0_8
8 real(dp), parameter :: epsrel = 1.e-3_8
9
10
11 real(dp), dimension(2) :: points
12 real(dp) :: Integ1
13
14 points(:2) = [1._8, sqrt(2._8)]
15
16 call qagp(fquad454, zero, 3._dp, points, integ1, epsabs, epsrel, err, neval=n)
17 print "(A)", 'integrate(x^3 \log{|(x^2-1)(x^2-2)|}, x, 0, 3) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
18 ! integrate(x^3 \log{|(x^2-1)(x^2-2)|}, x, 0, 3) = 52.7408058906746 (Error=0.0001754605945) with N = 777
19
20contains
21
23 function fquad454(x) result(y)
24 implicit none
25 real(dp) :: y
26 real(dp), intent(IN) :: x
27 y = x**3 * log(abs((x**2 - 1._dp) * (x**2 - 2._dp)))
28 end function fquad454
29
30end program ex_qagp

Subroutine QAWO for Oscillatory functions

Subroutine qawo computes the integral of a function f multiplied by a oscillatory function of $\cos{(\omega x)}$ (for flgw=1) or $\sin{(\omega x)}$ (for flgw=1),

Examples of use

1program ex_qawo
2 USE numfor, only: dp, zero, m_pi, str, qawo
3 implicit none
4 real(dp) :: err
5 integer :: N
6
7 real(dp) :: Integ1
8 ! Fifth argument, flgw=1 => weight function is cosine
9 ! fquad456 => log(x)
10 call qawo(fquad456, zero, 1._dp, 10 * m_pi, 2, integ1)
11 print "(A)", 'integrate(sin(10 pi x) log(x), x, 0, 1) = '//str(integ1)
12 ! integrate(sin(10 pi x) log(x), x, 0, 1) = -0.1281368491001
13
14
15 call qawo(fquad456, zero, 1._dp, 10 * m_pi, 1, integ1, epsabs=zero, epsrel=1.e-3_dp, abserr=err, neval=n)
16
17 print "(A)", 'integrate(cos(10 pi x) log(x), x, 0, 1) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
18 ! integrate(cos(10 pi x) log(x), x, 0, 1) = -0.0489888171135 (Error=4.8212107917056e-10) with N = 305
19contains
21 function fquad456(x) result(y)
22 implicit none
23 real(dp) :: y
24 real(dp), intent(IN) :: x
25 y = zero
26 IF (x > 0) y = log(x)
27 end function fquad456
28
29end program ex_qawo

Subroutine QAWF for Fourier Transforms

Subroutine qawf computes the (possibly semi-infinite) cosine (for flgw=1) or sine (for flgw=2) Fourier transform of f using an adaptive algorithm

Examples of use

1program ex_qawf
2 USE numfor, only: dp, zero, m_pi, str, qawf
3 implicit none
4 real(dp) :: err
5 integer :: N, ier
6
7 real(dp) :: Integ1
8 ! Fourth argument, flgw=1 => weight function is cosine
9 ! fquad457 => 1/\sqrt(x)
10 call qawf(fquad457, zero, m_pi / 2, 1, integ1)
11 print "(A)", 'integrate(cos(pi x /2 ) / sqrt(x), x, 0, inf ) = '//str(integ1)
12 ! integrate(cos(pi x /2 ) / sqrt(x), x, 0, inf ) = 0.9999999999321
13
14
15 call qawf(fquad457, zero, m_pi / 2, 1, integ1, epsabs=1.e-3_dp, abserr=err, neval=n, ier=ier)
16 print "(A)", 'integrate(cos(pi x /2 ) / sqrt(x), x, 0, inf ) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
17 ! integrate(cos(pi x /2 ) / sqrt(x), x, 0, inf ) = 0.9999969531148 (Error=0.0005923419178) with N = 380
18contains
20 function fquad457(x) result(y)
21 implicit none
22 real(dp) :: y
23 real(dp), intent(IN) :: x
24 y = 0._8
25 IF (x > 0) y = 1.0_8 / sqrt(x)
26 end function fquad457
27
28end program ex_qawf

Subroutine QAWS for integration of singular functions

Subroutine qaws computes the integral of a function f multiplied by a weight function with an algebraic-logarithmic singularity chosen according to the flag flgw,

flgw Weight function
1 $(x - a)^\alpha (b - x)^\beta$
2 $(x - a)^\alpha (b - x)^\beta \log{(x-a)}$
3 $(x - a)^\alpha (b - x)^\beta \log{(b-x)}$
4 $(x - a)^\alpha (b - x)^\beta \log{(x-a)} \log{(b-x)}$

The points $(a,b)$ are the limits of the integration domain.

Examples of use

1program ex_qaws
2 USE numfor, only: dp, zero, str, qaws
3 implicit none
4
5 real(dp) :: err
6 integer :: N
7
8 real(dp) :: Integ1
9 real(dp) :: alfa, beta
10 integer :: flgw
11 alfa = 0._dp
12 beta = 0._dp
13 flgw = 2 ! weight function => log(x-a)^alfa = log(x)
14 ! fquad458 = 1/(1 + \ln(x)^{2})^{2}
15 call qaws(fquad458, zero, 1._dp, alfa, beta, flgw, integ1)
16 print "(A)", 'integrate(log(x)/(1 + ln(x)^{2})^{2}, 0, 1) = '//str(integ1)
17 ! integrate(log(x)/(1 + ln(x)^{2})^{2}, 0, 1) = -0.1892750041577
18
19
20 ! Get estimation of error and number of evaluations
21 call qaws(fquad458, zero, 1._dp, alfa, beta, flgw, integ1, epsrel=1.e-3_dp, abserr=err, neval=n)
22 print "(A)", 'integrate(log(x)/(1 + ln(x)^{2})^{2}, 0, 1) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
23 ! integrate(log(x)/(1 + ln(x)^{2})^{2}, 0, 1) = -0.1892747444915 (Error=2.2226716553438e-06) with N = 40
24contains
25
27 function fquad458(x) result(y)
28 implicit none
29 real(dp) :: y
30 real(dp), intent(IN) :: x
31 y = 1 / (1 + log(x)**2)**2
32 end function fquad458
33
34end program ex_qaws

Subroutine QAWC for Cauchy principal values

Subroutine qawc computes Cauchy Principal Value of the integral of f using the bisection strategy of qag

Examples of use

1program ex_qaws
2 USE numfor, only: dp, str, qawc, nf_minf
3 implicit none
4 real(dp) :: err
5 integer :: N, ier
6
7 real(dp) :: Integ1
8 call qawc(fquad459, -1._dp, 5._dp, 0._dp, integ1)
9 print "(A)", 'integrate(1/(x(5 x^3 + 6), -1, 5)) = '//str(integ1)
10 ! integrate(1/(x(5 x^3 + 6), -1, 5)) = -0.0899440069576
11
12
13 call qawc(fquad459, -1._dp, 5._dp, 0._dp, integ1, epsrel=1.e-3_dp, abserr=err, neval=n, ier=ier)
14 print "(A)", 'integrate(1/(x(5 x^3 + 6), -1, 5)) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
15 ! integrate(1/(x(5 x^3 + 6), -1, 5)) = -0.0899440069584 (Error=1.1852922354142e-06) with N = 215
16contains
17
19 function fquad459(x) result(y)
20 implicit none
21 real(dp) :: y
22 real(dp), intent(IN) :: x
23 y = 1 / (5 * x**3 + 6)
24 end function fquad459
25
26end program ex_qaws

Use of the extra parameter info

Most adaptive routines allow for an optional input/output argument info. It is mainly used to pass information on details of the calculations. This argument is a variable of the type d_qp_extra for integration of real functions, and of the type c_qp_extra for integration of complex functions.

Following Piessend et al we illustrate one of its use considering the integral

\[ \int_{0}^{1} \left| x^{2} + 2 x -1 \right|^{-1/2} \, dx = \arcsin(1/\sqrt(3)) - \frac{\pi}{2} - \log(3)/2 \approx 1.504622 \]

In this case the integral converges

1program ex_qp_extra
2 USE numfor, only: dp, zero, m_pi_2, str, d_qp_extra, qag, qags
3 implicit none
4
5 real(dp) :: err
6 integer :: ier, N
7 integer :: j, k
8 type(d_qp_extra) :: info
9 real(dp) :: Integ1
10
11 ! Initialize allowing up to 50 subintervals
12 ! This is optional, by default it will allow for 500 subintervals
13 info = 50
14
15 call qags(fquad4510, zero, 1._dp, integ1, epsabs=zero, epsrel=1.e-4_dp, gkrule='qk15', abserr=err, neval=n, ier=ier, info=info)
16
17 print "(A)", 'integrate(1/sqrt(|x^2 + 2x -2|), x, 0, 1) = '//str(integ1)//" (Error="//str(err)//") with N = "//str(n)
18 print "(A)", "Difference = "//str(abs(integ1 - (m_pi_2 - asin(1 / sqrt(3._dp)) + log(3._dp) / 2)))//&
19 & " relative to analytical value"
20 print "(A)", "Error code = "//str(ier)
21 print "(A)", "Meaning:"
22 print "(A)", "-------"
23 print "(A)", info%msg
24 print "(A)", ""
25 print "(A)", " j Limits of subinterval S(j) Integral over S(j) Error on S(j) Error relative "
26 do k = 1, info%last
27 j = info%iord(k)
28 print "(I2, 5(es18.11,1x))", k, info%alist(j), info%blist(j), info%rlist(j), info%elist(j), info%elist(j) / info%rlist(j)
29 end do
30
31contains
32
34 function fquad4510(x) result(y)
35 implicit none
36 real(dp) :: y
37 real(dp), intent(IN) :: x
38 y = 1._dp / sqrt(abs(x**2 + 2 * x - 2))
39 end function fquad4510
40end program ex_qp_extra
41

giving as result information on the algorithm progress

integrate(1/sqrt(|x^2 + 2x -2|), x, 0, 1) = 1.5045599601491 (Error=0.0001175019052) with N = 735
Difference = 6.2802309494403e-05 relative to analytical value
Error code = 0
Meaning:
-------
Routine qags: Normal and reliable termination of the routine. It is assumed that the requested accuracy has been achieved.
j Limits of subinterval S(j) Integral over S(j) Error on S(j) Error relative
1 7.32050657272E-01 7.32050895691E-01 6.61788951586E-04 2.67199722355E-04 4.03753676628E-01
2 7.50000000000E-01 8.75000000000E-01 2.59695449021E-01 2.34942884257E-05 9.04686181996E-05
3 7.32055664062E-01 7.32116699219E-01 6.35457321524E-03 2.21241348822E-05 3.48160830520E-03
4 7.34375000000E-01 7.50000000000E-01 9.20419410896E-02 6.43439223180E-06 6.99071766157E-05
5 7.32050895691E-01 7.32051849365E-01 7.77806648660E-04 1.18010670163E-06 1.51722372605E-03
6 7.32421875000E-01 7.34375000000E-01 3.10999052060E-02 3.36607750891E-07 1.08234333405E-05
7 7.31445312500E-01 7.31933593750E-01 1.48085753181E-02 2.20888795975E-08 1.49162759570E-06
8 0.00000000000E+00 5.00000000000E-01 4.31717842526E-01 7.65314094754E-10 1.77271824179E-09
9 7.32051849365E-01 7.32055664062E-01 1.27128139071E-03 5.82935274600E-10 4.58541499043E-07
10 6.87500000000E-01 7.18750000000E-01 1.03290479678E-01 5.00387915491E-10 4.84447276312E-09
11 7.26562500000E-01 7.30468750000E-01 3.68841542687E-02 3.09541720462E-10 8.39226834936E-09
12 7.32040405273E-01 7.32048034668E-01 1.67639077543E-03 4.48780642525E-11 2.67706461466E-08
13 7.32177734375E-01 7.32421875000E-01 8.59296586961E-03 4.15971773602E-12 4.84084051903E-10
14 7.32048034668E-01 7.32049942017E-01 7.89650429625E-04 1.86687676455E-12 2.36418128137E-09
15 6.25000000000E-01 6.87500000000E-01 1.26121798856E-01 1.15514630686E-12 9.15897424031E-12
16 7.30468750000E-01 7.31445312500E-01 1.63018589630E-02 9.07839605128E-13 5.56893301058E-11
17 7.18750000000E-01 7.26562500000E-01 4.43801352397E-02 4.94821770138E-13 1.11496228541E-11
18 5.00000000000E-01 6.25000000000E-01 1.70177840734E-01 1.29207244330E-13 7.59248347332E-13
19 7.32025146484E-01 7.32040405273E-01 1.97766297002E-03 3.28124217391E-14 1.65915134361E-11
20 7.32050418854E-01 7.32050657272E-01 2.53372044356E-04 1.18649893348E-14 4.68283285354E-11
21 7.31994628906E-01 7.32025146484E-01 2.61075315192E-03 2.57842783013E-15 9.87618392123E-13
22 8.75000000000E-01 1.00000000000E+00 1.45769659141E-01 1.76224733434E-15 1.20892601707E-14
23 7.31933593750E-01 7.31994628906E-01 3.57974836972E-03 1.00512858960E-15 2.80781911405E-13
24 7.32049942017E-01 7.32050418854E-01 3.29764477394E-04 3.75908490859E-16 1.13993021271E-12
25 7.32116699219E-01 7.32177734375E-01 3.38357303442E-03 9.07748753959E-17 2.68281117247E-14

As we can see, the general integrator qags is able to integrate this singular function but needs several function evaluations, with an important error in the vicinity of $\sqrt{3}-1$. We can use this additional information and replace the routine qags by qagp. It only need to change line number XX to:

call qagp(fquad4510, Zero, 1._dp, [sqrt(3._dp) - 1], Integ1, epsabs=Zero, epsrel=1.e-4_dp, abserr=err, Neval=N, ier=ier, info=info)

obtaining a similar (and somewhat better) result with less subdivisions of the interval and function evaluations

integrate(1/sqrt(|x^2 + 2x -2|), x, 0, 1) = 1.5046227555966 (Error=3.5719196844752e-07) with N = 462
Difference = 6.8619860904562e-09 relative to analytical value
Error code = 0
Meaning:
-------
Routine qagp: Normal and reliable termination of the routine. It is assumed that the requested accuracy has been achieved.
j Limits of subinterval S(j) Integral over S(j) Error on S(j) Error relative
1 7.09174219832E-01 7.32050807569E-01 1.60069187812E-01 7.73744297312E-02 4.83381160290E-01
2 7.32050807569E-01 7.40424219832E-01 9.66937903941E-02 4.68487167016E-02 4.84505949251E-01
3 0.00000000000E+00 3.66025403784E-01 2.93171381140E-01 8.17539976481E-09 2.78860771915E-08
4 3.66025403784E-01 5.49038105677E-01 1.98297261193E-01 5.21098430467E-09 2.62786499083E-08
5 8.66025403784E-01 1.00000000000E+00 1.58478399902E-01 3.86710284589E-09 2.44014505969E-08
6 5.49038105677E-01 6.40544456623E-01 1.37342006395E-01 3.52212881075E-09 2.56449494455E-08
7 7.99038105677E-01 8.66025403784E-01 1.13597226144E-01 2.80948755739E-09 2.47320084544E-08
8 6.40544456623E-01 6.86297632096E-01 9.61465434121E-02 2.43826866463E-09 2.53599201604E-08
9 7.65544456623E-01 7.99038105677E-01 8.08861095193E-02 2.01475142213E-09 2.49084970721E-08
10 6.86297632096E-01 7.09174219832E-01 6.76512094521E-02 1.70646046002E-09 2.52243895393E-08
11 7.48797632096E-01 7.65544456623E-01 5.73966267734E-02 1.43490327440E-09 2.49997840476E-08
12 7.40424219832E-01 7.48797632096E-01 4.06573774170E-02 1.01831363346E-09 2.50462203455E-08