Numerical integration

Integration in scattering is used quite often. A good example is an averaging of a form-factor over certain parameter, like size of a particle, in small angle scattering. The best if a solution of this integral can be found analytically, but this is not always possible. Either integral doesn't have analytical solution or its solution is too specific and can be applied for a particular limited cases only, like particular distribution function.

Also, analytical calculation could require a lot of time and result usually is needed as soon as possible. Plus, one can use a numerical integration for a proof of a found analytical solution.

Here, methods of numerical integration of a functor object are presented. Currently, numerical integration routines in ESCAPE are based on Gauss-Kronrod Quadratures method. Below are few examples.

Standard integration

ESCAPE math module supports two types of integration with differential of functor variable or functor parameter.

To find a numerical solution of integral of type $\int_a^bf(x, p)dx$ or $\int_a^bf(x, p)dp$, where $x$ is a variable of functor $f(x, p)$ and $p$ is a parameter and $a, b$ are the finite limits of integration, one should use integral method as following:

integral(F, vp, a, b, storage_size=7, epsabs=None,
              epsrel=None, maxiter=None)

where $F$ - functor with double return type, $x$ - integration variable or parameter, $a, b$ - integration limits, storage_size - number of quadrature points [7, 15, 21, 51, 61]), epsabs and epsrel - absolute and relative tolerances for adaptive integration, maxiter - maximum number of iterations of adaptive integration. if maxiter is zero, adaptive integration is switched off. If you integrate a functor of one variable $x$ over this variable, the result if integrate method is a dependent parameter. if you integrate a functor of two variables over one of them like $\int F(x, y)dx$, the returned object will be a functor $I(y)$. Below there are several examples.

Averaging of a function with respect to distribution function

ESCAPE implements several useful well-known distribution functions and their integrals. For the integration we still use Gauss-Kronrod quadratures and there is no support for infinite integration limits, but for typical distribution function we present here, the integration limits either are finite or can be replaced by finite integration limits still preserving high accuracy of the integration result. As for the standard integration routine, there is a support for integration over a functor variable or over a parameter.

Integration over variable

Integration over variable can be described by the following equation:

$F(x_0)=\int_a^b f(x)G(x, x_0; \sigma(x_0))dx$, where $x_0$ - is a variable which correspond to the mean value, $G(x, x_0; \sigma(x_0))$ - distribution function of two variables $x$ and $x_0$, $\sigma(x_0)$ - $\sigma$ or FWHM functor, $a$ and $b$ are the integration limits, which depend on distribution function and its standard deviation.

Integration over parameter

Integration over parameter can be described by the following equation:

$F(x; p_0)=\int_a^b f(x; p)G(x, p_0; \sigma(p_0))dp$, where $p_0$ - is a parameter which correspond to the mean value, $p$ - integration parameter, $G(p, p_0; \sigma(p_0))$ - distribution functor of one variable and $p_0$ parameter as a mean value, $\sigma(p_0)$ - $\sigma$ or FWHM parameter.

As an example we integrate a function $f(x)=sin^2(px)$ over the variable $x$ and over the parameter $p$. The parameter $p$ defines the oscillation frequency of the integrated functor. At high oscillation frequencies, non-adaptive Gauss-Kronrod method can fail. That's why in all examples below we have set the adaptive integration method with quite high maxiter=100 parameter, just be sure that adaptive integration has large enough maximum number of iterations and will converge to the proper solution.

Gamma distribution function

We use the following form of the gamma distribution function:

$G(x; x_0, \sigma)=1/\theta\exp((k-1)\log(t)-x/\theta-\log(\Gamma(k)))$

where

$\theta = \sigma^2 x_0$

$k = x_0 / \theta$

$t = x / \theta$

Schulz-Zimm distribution function

$G(x; x_0, \sigma)=t / x \exp((k - 1) \log(t) - t - \log(\Gamma(k)))$

where

$k = 1/\sigma^2$

$t = kx/x_0$

LogNorm distribution function

$G(x; x_0, \sigma)=\exp(\log^2(x / x_0) / 2 / \sigma^2))/(\sqrt{2\pi}\sigma x)$

Normal distribution function

$G(x; x_0, \sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{(-(x - x_0)^2 / (2 \sigma^2))}$

Triangular distribution function

$G(x; x_0, \sigma) = 2 (x - a) / (b - a) / (x_0 - a) $,

for $a<=x<=x_0$

$G(x; x_0, \sigma) = 2 (b - x) / (b - a) / (b - x_0) $,

for $x_0<=x<=b$,

where $a=x_0-\sigma$ and $b=x_0+\sigma$

Uniform distribution function

if $w$ is a FWHM.

$G(x; x_0, \sigma) = 1.0 / w $,

for $x_0-w/2<=x<=x_0+w/2$

Custom distribution function

Average method allows one to find an average of any functor with user-defined distribution function. This is how the code for the user-defined Gaussian distribution looks like.