SF_INTEGRATE

Description

SciFortran module for function integration

Quick access

Routines:: get_quadrature_weights(), kronig(), quad(), simps(), simps2d(), trapz(), trapz2d()

Used modules

gauss_quadrature

Subroutines and functions

interface sf_integrate/trapz(f, a, b, dh, x[, a, b, n])

This function evaluates the integral of a discretized real or complex function using the trapezoidal rule (2nd order), or of a function passed as argument. The specific procedures cover

a sampled function on a uniform grid defined by the limits a, b
a sampled function on a uniform grid defined by the step dh
a sampled function on a non-uniform grid
a function evaluated between the limits a, b
a function evaluated on a non-uniform grid

The quadrature weights are provided by get_quadrature_weights with nrk=2.

Parameters:

f (various shapes) [real/complex] – integrand: function f(x), or array of its samples f(x_i)
a [real]
b [real] – integration limits a (lower), b (upper); default 0,1 for functions
dh [real] – uniform grid spacing
x (various shapes) [real]

Result:

sum [real/complex]
int [real/complex]

Options:

n [integer] – number of grid points (default 100)

interface sf_integrate/trapz2d(func, xrange, yrange, nx, ny[, n0, iterative, threshold, dhx, dhy, xrange, yrange])

This function evaluates the 2-dimensional integral of a real or complex function using the trapezoidal rule (2nd order). The specific procedures cover a function of two variables, a function evaluated recursively as a nested 1-dimensional integral, and a sampled function on a 2-dimensional grid.

Parameters:

func (various shapes) [real/complex, external] – integrand: function f(x) with x=[x,y], or 2d array of samples f(x_i,y_j)
xrange (2) [real]
yrange (2) [real] – integration limits [lower,upper] along x and y (default [0,1] for samples)
nx [integer]
ny [integer] – number of grid points along x and y

Result:

int [real/complex]

Options:

n0 [integer] – initial number of grid points along x and y (default 51)
iterative [logical] – if T refine grid N->2*N-10 until relative change < threshold (default F)
threshold [real] – relative accuracy of the grid refinement (default 5d-3 if iterative)
dhx [real]
dhy [real] – grid spacings along x and y (default: from xrange,yrange and grid size)

interface sf_integrate/simps(f, dh, a, b, x[, a, b, n])

This function evaluates the integral of a discretized real or complex function using Simpson’s rule (4th order), or of a function passed as argument. The specific procedures cover

a sampled function on a uniform grid defined by the step dh
a sampled function on a uniform grid defined by the limits a, b
a sampled function on a non-uniform grid
a function evaluated between the limits a, b
a function evaluated on a non-uniform grid

The quadrature weights are provided by get_quadrature_weights with nrk=4.

Parameters:

f (various shapes) [real/complex] – integrand: function f(x), or array of its samples f(x_i)
dh [real] – uniform grid spacing
a [real]
b [real] – integration limits a (lower), b (upper); default 0,1 for functions
x (various shapes) [real]

Result:

sum [real/complex]
int [real/complex]

Options:

n [integer] – number of grid points (default 100)

interface sf_integrate/simps2d(func, xrange, yrange, nx, ny[, n0, iterative, threshold, dhx, dhy, xrange, yrange])

This function evaluates the 2-dimensional integral of a real or complex function using Simpson’s rule (4th order). The specific procedures cover a function of two variables, a function evaluated recursively as a nested 1-dimensional integral, and a sampled function on a 2-dimensional grid.

Parameters:

func (various shapes) [real/complex, external] – integrand: function f(x) with x=[x,y], or 2d array of samples f(x_i,y_j)
xrange (2) [real]
yrange (2) [real] – integration limits [lower,upper] along x and y (default [0,1] for samples)
nx [integer]
ny [integer] – number of Simpson panels along x and y (2*Nx, 2*Ny grid points)

Result:

int [real/complex]

Options:

n0 [integer] – initial number of Simpson panels along x and y (default 51)
iterative [logical] – if T refine grid N->2*N-10 until relative change < threshold (default F)
threshold [real] – relative accuracy of the grid refinement (default 1d-3 if iterative)
dhx [real]
dhy [real] – grid spacings along x and y (default: from xrange,yrange and grid size)

interface sf_integrate/quad(func, a, result, fsample, b[, b, epsabs, epsrel, key, inf, singular_endpoint, singular_points, cpole, alfa, beta, omega, weight_func, verbose, strict, ninterp])

This function evaluates the 1-dimensional integral of a function, or of a sampled function, using the adaptive routines of the QUADPACK library. The routine is chosen according to the type of integrand and integration region, as described in the guide at the top of the module (for example QAGS, QAGI, QAWO, QAWF, QAWC, QAWS). The optional argument weight_func selects the weight function \(W(x)\) of the integrand \(W(x) f(x)\) in the routines that integrate with a weight.

For QAWO and QAWF, which also require omega, the accepted values are

weight_func=1: \(W(x) = \cos(\omega x)\)
weight_func=2: \(W(x) = \sin(\omega x)\)

For QAWS, which also requires alfa and beta, the accepted values are

weight_func=1: \(W(x) = (x-a)^{\alpha} (b-x)^{\beta}\)
weight_func=2: \(W(x) = (x-a)^{\alpha} (b-x)^{\beta} \ln(x-a)\)
weight_func=3: \(W(x) = (x-a)^{\alpha} (b-x)^{\beta} \ln(b-x)\)
weight_func=4: \(W(x) = (x-a)^{\alpha} (b-x)^{\beta} \ln(x-a) \ln(b-x)\)

with \(\alpha, \beta > -1\). The program stops if weight_func is outside these ranges. It is ignored by the other routines.

Parameters:

func [external] – integrand f(x)
a [real] – lower integration limit (finite limit of the range if b is absent)
result [real]
fsample (•) [real] – integrand sampled on a uniform grid of size(fsample) points in [a,b]

Options:

b [real] – upper integration limit (absent: semi-infinite range, needs inf or omega)
epsabs [real] – absolute accuracy requested (default 1d-12)
epsrel [real] – relative accuracy requested (default 1d-6)
key [integer] – Gauss-Kronrod rule of QAG: 1..6 = 15,21,31,41,51,61 points
inf [integer] – range of QAGI: -1 (-inf:a), 1 (a:inf), 2 (-inf:inf)
singular_endpoint [logical] – T if singular_endpoint exists (QAGS)
singular_points (•) [real] – location of singular points in QAGP
cpole [real] – location of the pole in QAWC f(x)/x-cpole
alfa [real]
beta [real] – exponents of the QAWS weight (x-a)**alfa*(b-x)**beta, both > -1
omega [real] – frequency of the weight function cos/sin(omega*x) in QAWO,QAWF
weight_func [integer] – weight function index
verbose [logical] – if T print the selected routine and the integration report (default F)
strict [logical] – T (default) stop on QUADPACK error code IER (only if verbose=F)
ninterp [integer] – order of the polynomial interpolation of fsample (default 3)

function sf_integrate/kronig(fi, wr, m)

This function evaluates the Kramers-Kronig transform of a function of a real variable. Given the values fi of the imaginary part on the uniform grid wr of M points, it returns the real part fr on the same grid,

\(f_r(w_i) = \frac{1}{\pi} \, \mathcal{P} \int_{w_1}^{w_M} \frac{f_i(w')}{w'-w_i} \, dw'\)

The principal value is evaluated by subtracting the singularity,

\(\mathcal{P} \int \frac{f_i(w')}{w'-w_i} dw' = \int \frac{f_i(w')-f_i(w_i)}{w'-w_i} dw' + f_i(w_i) \, \mathcal{P} \int \frac{dw'}{w'-w_i}\)

where the first integral is a sum with weight dh=wr(2)-wr(1), using the finite-difference derivative of \(f_i\) at the point \(w'=w_i\), and the second is known analytically, \(\ln\left(\frac{w_M-w_i}{w_i-w_1}\right)\). This term is set to zero at the end points \(i=1,M\), where it diverges. The cost scales as \(M^2\).

Parameters:

fi (m) [real]
wr (m) [real]

Options:

m [integer]

Result:

fr (m) [real]

subroutine sf_integrate/get_quadrature_weights(wt[, nrk])

This subroutine returns the quadrature weights wt for the integration of a function sampled on N=size(wt) points with unit step, such that

\(\int f(x) dx \approx dh \sum_{i=1}^{N} wt(i) \cdot f(x_i)\)

The order of the rule is selected by the optional argument nrk

nrk=2: trapezoidal rule (2nd order), \(wt = [1/2, 1, ..., 1, 1/2]\)
nrk=4: Simpson’s rule (4th order), the default

For nrk=4 the weights depend on N

N=1: \(wt = 1\)
N=2: \(wt = [1/2, 1/2]\)
N=3: Simpson’s rule \([1/3, 4/3, 1/3]\)
N=4: Simpson’s 3/8 rule \([3/8, 9/8, 9/8, 3/8]\)
N>=5 odd: composite Simpson’s rule \([1/3, 4/3, 2/3, ..., 4/3, 1/3]\)
N>=6 even: composite Simpson’s rule, closed by a Simpson’s 3/8 rule on the last four points

The program stops if nrk is neither 2 nor 4.

Parameters:: wt (•) [real] – Quadrature weights
Options:: nrk [integer] – Order of the rule