quadrature_module Module

  1. quadrature_module.F90
  2. quadrature_module

Integration of functions using adaptive Guassian quadrature.

Author

  • Jacob Williams


Uses

  • module~~quadrature_module~~UsesGraph module~quadrature_module quadrature_module iso_fortran_env iso_fortran_env module~quadrature_module->iso_fortran_env

Contents

Variables

Abstract Interfaces

Derived Types

Functions

Subroutines


Variables

TypeVisibilityAttributesNameInitial
integer, public, parameter:: quadrature_wp =real64

real kind used by this module [8 bytes]
integer, private, parameter:: wp =quadrature_wp

local copy of fmin_rk with a shorter name
real(kind=wp), private, parameter:: zero =0.0_wp
real(kind=wp), private, parameter:: one_half =0.5_wp
real(kind=wp), private, parameter:: one =1.0_wp
real(kind=wp), private, parameter:: two =2.0_wp
real(kind=wp), private, parameter:: three =3.0_wp
real(kind=wp), private, parameter:: four =4.0_wp
type(quadrature_method), private, parameter:: quad_gauss_6 =quadrature_method(6, 'Adaptive 6-point Legendre-Gauss')
type(quadrature_method), private, parameter:: quad_gauss_8 =quadrature_method(8, 'Adaptive 8-point Legendre-Gauss')
type(quadrature_method), private, parameter:: quad_gauss_10 =quadrature_method(10, 'Adaptive 10-point Legendre-Gauss')
type(quadrature_method), private, parameter:: quad_gauss_12 =quadrature_method(12, 'Adaptive 12-point Legendre-Gauss')
type(quadrature_method), private, parameter:: quad_gauss_14 =quadrature_method(14, 'Adaptive 14-point Legendre-Gauss')
type(quadrature_method), public, parameter, dimension(5):: set_of_quadrature_methods =[quad_gauss_6, quad_gauss_8, quad_gauss_10, quad_gauss_12, quad_gauss_14]

Abstract Interfaces

abstract interface

  • private function func_1d(me, x) result(f)

    1d user function f(x)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_1d), intent(inout) :: me
    real(kind=wp), intent(in) :: x

    Return Value real(kind=wp)

abstract interface

  • private function func_2d(me, x, y) result(f)

    2d user function f(x,y)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_2d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: y

    Return Value real(kind=wp)

abstract interface

  • private function func_3d(me, x, y, z) result(f)

    3d user function f(x,y,z)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_3d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: y
    real(kind=wp), intent(in) :: z

    Return Value real(kind=wp)

abstract interface

  • private function func_4d(me, x, y, z, q) result(f)

    4d user function f(x,y,z,q)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_4d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: y
    real(kind=wp), intent(in) :: z
    real(kind=wp), intent(in) :: q

    Return Value real(kind=wp)

abstract interface

  • private function func_5d(me, x, y, z, q, r) result(f)

    5d user function f(x,y,z,q,r)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_5d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: y
    real(kind=wp), intent(in) :: z
    real(kind=wp), intent(in) :: q
    real(kind=wp), intent(in) :: r

    Return Value real(kind=wp)

abstract interface

  • private function func_6d(me, x, y, z, q, r, s) result(f)

    6d user function f(x,y,z,q,r,s)

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_6d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: y
    real(kind=wp), intent(in) :: z
    real(kind=wp), intent(in) :: q
    real(kind=wp), intent(in) :: r
    real(kind=wp), intent(in) :: s

    Return Value real(kind=wp)

abstract interface

  • private function gauss_func(me, x, h) result(f)

    guass quadrature formula

    Arguments

    TypeIntentOptionalAttributesName
    class(integration_class_1d), intent(inout) :: me
    real(kind=wp), intent(in) :: x
    real(kind=wp), intent(in) :: h

    Return Value real(kind=wp)


Derived Types

type, public :: quadrature_method

quadrature methods

Components

TypeVisibilityAttributesNameInitial
integer, public :: n_points =0
character(len=100), public :: name =''

type, private, abstract :: integration_class

type, public, extends(integration_class) :: integration_class_1d

single integration class: for 1d integration of f(x)

Components

TypeVisibilityAttributesNameInitial
procedure(func_1d), public, pointer:: fun=> null()

function f(x) to be integrated
procedure(gauss_func), public, pointer:: g=> null()

the guass quadrature formula to use
real(kind=wp), public :: a =zero

lower limit of integration
real(kind=wp), public :: b =zero

upper limit of integration (may be less than a)
real(kind=wp), public :: tol =zero

the requested relative error tolerance.
real(kind=wp), public :: val =zero

the value of x. Only used for multiple integration to pass to the inner integrals

Type-Bound Procedures

procedure, public :: dgauss_generic

core integration routine. refactored from SLATEC with selectable quadrature method

procedure, public :: initialize => initialize_integration_class

to set up the class

procedure, public :: integrate => integrate_1d

to integrate the function fun

type, public, extends(integration_class) :: integration_class_2d

double integration class: for 2d integration of f(x,y)

Components

TypeVisibilityAttributesNameInitial
procedure(func_2d), public, pointer:: fxy=> null()

function f(x,y) to be integrated
type(integration_class_1d), public :: iy

for the dy integration
type(integration_class_1d), public :: ix

for the dx integration

Type-Bound Procedures

procedure, public :: initialize => initialize_integration_class_2d

to set up the class

procedure, public :: integrate => integrate_2d

to integrate the function fxy

type, public, extends(integration_class) :: integration_class_3d

double integration class: for 3d integration of f(x,y,z)

Components

TypeVisibilityAttributesNameInitial
procedure(func_3d), public, pointer:: fxyz=> null()

function f(x,y,z) to be integrated
type(integration_class_1d), public :: iz

for the dz integration
type(integration_class_1d), public :: iy

for the dy integration
type(integration_class_1d), public :: ix

for the dx integration

Type-Bound Procedures

procedure, public :: initialize => initialize_integration_class_3d

to set up the class

procedure, public :: integrate => integrate_3d

to integrate the function fxyz

type, public, extends(integration_class) :: integration_class_4d

double integration class: for 4d integration of f(x,y,z,q)

Components

TypeVisibilityAttributesNameInitial
procedure(func_4d), public, pointer:: fxyzq=> null()

function f(x,y,z,q) to be integrated
type(integration_class_1d), public :: iq

for the dq integration
type(integration_class_1d), public :: iz

for the dz integration
type(integration_class_1d), public :: iy

for the dy integration
type(integration_class_1d), public :: ix

for the dx integration

Type-Bound Procedures

procedure, public :: initialize => initialize_integration_class_4d

to set up the class

procedure, public :: integrate => integrate_4d

to integrate the function fxyzq

type, public, extends(integration_class) :: integration_class_5d

double integration class: for 5d integration of f(x,y,z,q,r)

Components

TypeVisibilityAttributesNameInitial
procedure(func_5d), public, pointer:: fxyzqr=> null()

function f(x,y,z,q,r) to be integrated
type(integration_class_1d), public :: ir

for the dr integration
type(integration_class_1d), public :: iq

for the dq integration
type(integration_class_1d), public :: iz

for the dz integration
type(integration_class_1d), public :: iy

for the dy integration
type(integration_class_1d), public :: ix

for the dx integration

Type-Bound Procedures

procedure, public :: initialize => initialize_integration_class_5d

to set up the class

procedure, public :: integrate => integrate_5d

to integrate the function fxyzqr

type, public, extends(integration_class) :: integration_class_6d

double integration class: for 6d integration of f(x,y,z,q,r,s)

Components

TypeVisibilityAttributesNameInitial
procedure(func_6d), public, pointer:: fxyzqrs=> null()

function f(x,y,z,q,r,s) to be integrated
type(integration_class_1d), public :: is

for the ds integration
type(integration_class_1d), public :: ir

for the dr integration
type(integration_class_1d), public :: iq

for the dq integration
type(integration_class_1d), public :: iz

for the dz integration
type(integration_class_1d), public :: iy

for the dy integration
type(integration_class_1d), public :: ix

for the dx integration

Type-Bound Procedures

procedure, public :: initialize => initialize_integration_class_6d

to set up the class

procedure, public :: integrate => integrate_6d

to integrate the function fxyzqrs

Functions

private function g6(me, x, h) result(f)

6-point method.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: h

Return Value real(kind=wp)

private function g8(me, x, h) result(f)

This is the 8-point formula from the original SLATEC routine DGAUS8.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: h

Return Value real(kind=wp)

private function g10(me, x, h) result(f)

10-point method.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: h

Return Value real(kind=wp)

private function g12(me, x, h) result(f)

12-point method.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: h

Return Value real(kind=wp)

private function g14(me, x, h) result(f)

14-point method.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in) :: h

Return Value real(kind=wp)


Subroutines

private subroutine initialize_integration_class(me, fx, xl, xu, tolx, methodx)

Initialize the 1D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
procedure(func_1d) :: fx

1d function: f(x)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
integer, intent(in) :: methodx

quadrature method to use for x

private subroutine initialize_integration_class_2d(me, fxy, xl, xu, yl, yu, tolx, toly, methodx, methody)

Initialize the 2D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_2d), intent(inout) :: me
procedure(func_2d) :: fxy

2d function: f(x,y)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: yl

y integration lower bound
real(kind=wp), intent(in) :: yu

y integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
real(kind=wp), intent(in) :: toly

error tolerance for dy integration
integer, intent(in) :: methodx

quadrature method to use for x
integer, intent(in) :: methody

quadrature method to use for y

private subroutine initialize_integration_class_3d(me, fxyz, xl, xu, yl, yu, zl, zu, tolx, toly, tolz, methodx, methody, methodz)

Initialize the 3D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_3d), intent(inout) :: me
procedure(func_3d) :: fxyz

3d function: f(x,y,z)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: yl

y integration lower bound
real(kind=wp), intent(in) :: yu

y integration upper bound
real(kind=wp), intent(in) :: zl

z integration lower bound
real(kind=wp), intent(in) :: zu

z integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
real(kind=wp), intent(in) :: toly

error tolerance for dy integration
real(kind=wp), intent(in) :: tolz

error tolerance for dz integration
integer, intent(in) :: methodx

quadrature method to use for x
integer, intent(in) :: methody

quadrature method to use for y
integer, intent(in) :: methodz

quadrature method to use for z

private subroutine initialize_integration_class_4d(me, fxyzq, xl, xu, yl, yu, zl, zu, ql, qu, tolx, toly, tolz, tolq, methodx, methody, methodz, methodq)

Initialize the 4D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_4d), intent(inout) :: me
procedure(func_4d) :: fxyzq

4d function: f(x,y,z,q)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: yl

y integration lower bound
real(kind=wp), intent(in) :: yu

y integration upper bound
real(kind=wp), intent(in) :: zl

z integration lower bound
real(kind=wp), intent(in) :: zu

z integration upper bound
real(kind=wp), intent(in) :: ql

q integration lower bound
real(kind=wp), intent(in) :: qu

q integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
real(kind=wp), intent(in) :: toly

error tolerance for dy integration
real(kind=wp), intent(in) :: tolz

error tolerance for dz integration
real(kind=wp), intent(in) :: tolq

error tolerance for dq integration
integer, intent(in) :: methodx

quadrature method to use for x
integer, intent(in) :: methody

quadrature method to use for y
integer, intent(in) :: methodz

quadrature method to use for z
integer, intent(in) :: methodq

quadrature method to use for q

private subroutine initialize_integration_class_5d(me, fxyzqr, xl, xu, yl, yu, zl, zu, ql, qu, rl, ru, tolx, toly, tolz, tolq, tolr, methodx, methody, methodz, methodq, methodr)

Initialize the 5D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_5d), intent(inout) :: me
procedure(func_5d) :: fxyzqr

5d function: f(x,y,z,q,r)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: yl

y integration lower bound
real(kind=wp), intent(in) :: yu

y integration upper bound
real(kind=wp), intent(in) :: zl

z integration lower bound
real(kind=wp), intent(in) :: zu

z integration upper bound
real(kind=wp), intent(in) :: ql

q integration lower bound
real(kind=wp), intent(in) :: qu

q integration upper bound
real(kind=wp), intent(in) :: rl

r integration lower bound
real(kind=wp), intent(in) :: ru

r integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
real(kind=wp), intent(in) :: toly

error tolerance for dy integration
real(kind=wp), intent(in) :: tolz

error tolerance for dz integration
real(kind=wp), intent(in) :: tolq

error tolerance for dq integration
real(kind=wp), intent(in) :: tolr

error tolerance for dr integration
integer, intent(in) :: methodx

quadrature method to use for x
integer, intent(in) :: methody

quadrature method to use for y
integer, intent(in) :: methodz

quadrature method to use for z
integer, intent(in) :: methodq

quadrature method to use for q
integer, intent(in) :: methodr

quadrature method to use for r

private subroutine initialize_integration_class_6d(me, fxyzqrs, xl, xu, yl, yu, zl, zu, ql, qu, rl, ru, sl, su, tolx, toly, tolz, tolq, tolr, tols, methodx, methody, methodz, methodq, methodr, methods)

Initialize the 6D integration class. Must be called before integration is performed.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_6d), intent(inout) :: me
procedure(func_6d) :: fxyzqrs

6d function: f(x,y,z,q,r,s)
real(kind=wp), intent(in) :: xl

x integration lower bound
real(kind=wp), intent(in) :: xu

x integration upper bound
real(kind=wp), intent(in) :: yl

y integration lower bound
real(kind=wp), intent(in) :: yu

y integration upper bound
real(kind=wp), intent(in) :: zl

z integration lower bound
real(kind=wp), intent(in) :: zu

z integration upper bound
real(kind=wp), intent(in) :: ql

q integration lower bound
real(kind=wp), intent(in) :: qu

q integration upper bound
real(kind=wp), intent(in) :: rl

r integration lower bound
real(kind=wp), intent(in) :: ru

r integration upper bound
real(kind=wp), intent(in) :: sl

s integration lower bound
real(kind=wp), intent(in) :: su

s integration upper bound
real(kind=wp), intent(in) :: tolx

error tolerance for dx integration
real(kind=wp), intent(in) :: toly

error tolerance for dy integration
real(kind=wp), intent(in) :: tolz

error tolerance for dz integration
real(kind=wp), intent(in) :: tolq

error tolerance for dq integration
real(kind=wp), intent(in) :: tolr

error tolerance for dr integration
real(kind=wp), intent(in) :: tols

error tolerance for ds integration
integer, intent(in) :: methodx

quadrature method to use for x
integer, intent(in) :: methody

quadrature method to use for y
integer, intent(in) :: methodz

quadrature method to use for z
integer, intent(in) :: methodq

quadrature method to use for q
integer, intent(in) :: methodr

quadrature method to use for r
integer, intent(in) :: methods

quadrature method to use for s

private subroutine integrate_1d(me, ans, ierr, err)

Perform the 1D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private subroutine integrate_2d(me, ans, ierr, err)

Perform the 2D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_2d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private subroutine integrate_3d(me, ans, ierr, err)

Perform the 3D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_3d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private subroutine integrate_4d(me, ans, ierr, err)

Perform the 4D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_4d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private subroutine integrate_5d(me, ans, ierr, err)

Perform the 5D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_5d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private subroutine integrate_6d(me, ans, ierr, err)

Perform the 6D integration.

Arguments

TypeIntentOptionalAttributesName
class(integration_class_6d), intent(inout) :: me
real(kind=wp), intent(out) :: ans
integer, intent(out) :: ierr
real(kind=wp), intent(out) :: err

private recursive subroutine dgauss_generic(me, lb, ub, error_tol, ans, ierr, err)

Integrate a real function of one variable over a finite interval using the specified adaptive algorithm. Intended primarily for high accuracy integration or integration of smooth functions.

Read more…

Arguments

TypeIntentOptionalAttributesName
class(integration_class_1d), intent(inout) :: me
real(kind=wp), intent(in) :: lb

lower bound of the integration
real(kind=wp), intent(in) :: ub

upper bound of the integration
real(kind=wp), intent(in) :: error_tol

is a requested pseudorelative error tolerance. normally pick a value of abs(error_tol) so that dtol < abs(error_tol) <= 1.0e-3 where dtol is the larger of 1.0e-18 and the real(wp) unit roundoff d1mach(4). ans will normally have no more error than abs(error_tol) times the integral of the absolute value of fun(x). usually, smaller values of error_tol yield more accuracy and require more function evaluations.
real(kind=wp), intent(out) :: ans

computed value of integral
integer, intent(out) :: ierr

status code:

  • normal codes:
  • 1 : ans most likely meets requested error tolerance, or lb=ub.
  • -1 : lb and ub are too nearly equal to allow normal integration. ans is set to zero.
  • abnormal code:
  • 2 : ans probably does not meet requested error tolerance.
real(kind=wp), intent(out) :: err

an estimate of the absolute error in ans. the estimated error is solely for information to the user and should not be used as a correction to the computed integral.