In [1]:
import escape as esc
from escape.utils.widgets import show
import numpy as np
esc.require('0.9.7')

Loading material database from /home/dkor/Data/Development/workspace_escape/escape/python/src/escape/scattering/../data/mdb/materials.db

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, numpoints=7, epsabs=None,
              epsrel=None, maxiter=None)

where $F$ - functor with double return type, $x$ - integration variable or parameter, $a, b$ - integration limits, numpoints - 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.

In [2]:
X=esc.var("X")
Y=esc.var("Y")
In [3]:
p = esc.par("P", 1.0)
F1 = esc.sin(p*X)
In [4]:
#we create an integral of f(x; p) - I, which is a parameter  
I = esc.integral(F1, X, 0, 1.0)
print (type(I))

<class 'escape.core.objects.parameter_obj'>
In [5]:
show(I)
In [6]:
#integration of 2d functor
F2 = esc.sin(p*X)*esc.cos(p*Y)
I=esc.integral(F2, X, 0.0, 1.0, epsabs=1e-4, epsrel=1e-5, maxiter=10)
show(I, xlabel="Y", ylabel="I")
In [7]:
#integration over parameter
p = esc.par("P", 1.0)
F3 = esc.sin(p*X)
I = esc.integral(F3, p, 0, 1.0)
show(I)

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.

In [8]:
X0=esc.var("X0")

#Parameter 'p' and its mean value 'p0'
p=esc.par("p", 5.0, userlim=[0, 10])
#p0 is used only for integration over parameter
p0=esc.par("p0", 5.0, userlim=[0, 10])
#Sigma parameter used for gamma and schulz distribution functions
sigma=esc.par("Sigma", 0.2)

#functor to integrate
F=esc.pow(esc.sin(X*p), 2.0)

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$

In [9]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.gamma("G", X, p0, sigma)
show(G, coordinates=x, title="Gamma distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 10 iterations and 7 standard deviations for the limits
I = esc.average_gamma(F, sigma, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_gamma(F, sigma, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

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$

In [10]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.schulz("G", X, p0, sigma)
show(G, coordinates=x, title="Schulz-Zimm distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 100 iterations and 7 standard deviations for the limits
I = esc.average_schulz(F, sigma, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_schulz(F, sigma, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

LogNorm distribution function¶

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

In [11]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.lognorm("G", X, p0, sigma)
show(G, coordinates=x, title="LogNorm distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 10 iterations and 7 standard deviations for the limits
I = esc.average_lognorm(F, sigma, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_lognorm(F, sigma, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

Normal distribution function¶

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

In [12]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.normal("G", X, p0, sigma)
show(G, coordinates=x, title="Normal distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 10 iterations and 7 standard deviations for the limits
I = esc.average_normal(F, sigma*2.355, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_normal(F, sigma*2.355, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

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$

In [13]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.triangular("G", X, p0, sigma)
show(G, coordinates=x, title="Triangular distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 10 iterations and 7 standard deviations for the limits
I = esc.average_triangular(F, sigma*2, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_triangular(F, sigma*2, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

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$

In [14]:
# coordinates array
x=np.arange(0, 10, 0.01)

G=esc.uniform("G", X, p0, sigma)
show(G, coordinates=x, title="Uniform distribution")

#integration over variable
show(F, coordinates=x, title="F(x)=sin(px)^2")

#we use adaptive integration with max 100 iterations 
I = esc.average_uniform(F, sigma*2, X, X0, maxiter=100, epsrel=1e-8, epsabs=1e-8)
show(I, coordinates=x, title="Integration over variable", xlabel="X0", ylabel="I")

#integration over parameter
I = esc.average_uniform(F, sigma*2, p, p0, maxiter=100, epsrel=1e-8, epsabs=1e-8)
show(I, coordinates=x, title="Integration over parameter", xlabel="X", ylabel="I")

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.

In [15]:
# averaging over variable
G=1/np.sqrt(2*np.pi)/sigma*esc.exp(-(X-X0)*(X-X0)/(2*sigma*sigma))
G.variables
Out[15]:
[variable(name='X'), variable(name='X0')]
In [16]:
I=esc.average(F, G, X, X0, X0-10*sigma, X0+10*sigma, maxiter=100)
In [17]:
show(I, coordinates=np.arange(0, 10, 0.01), xlabel="X0")
In [18]:
# averaging over parameter
G=1/np.sqrt(2*np.pi)/sigma*esc.exp(-(p-p0)*(p-p0)/(2*sigma*sigma))
In [19]:
I=esc.average(F, G, p, p0, p0-5*sigma, p0+5*sigma, maxiter=100)
show(I, ylabel="I")
In [ ]: