import escape as esc
import numpy as np
esc.require("0.9.8")
Loading material database from C:\dev\escape-core\python\src\escape\scattering\..\data\mdb\materials.db
Numerical integration¶
Integration is common in scattering: for example, averaging a form factor over a size distribution in small-angle scattering, or over an orientation distribution. Physically, this is how ESCAPE accounts for sample polydispersity or for a quantity that is not fixed but distributed in the experiment. When no closed-form integral exists or the analytical solution is too narrow, numerical integration is used; it can also validate an analytical result.
ESCAPE provides numerical integration of functors over a variable or over a parameter (e.g. with a distribution), using Gauss–Kronrod quadrature with optional adaptive refinement (tolerances epsabs, epsrel, maxiter). Return type: integrating over the only variable gives a dependent parameter; integrating over one variable of a multi-variable functor gives a functor of the remaining variables; integrating over a parameter gives a functor or parameter depending on what variables remain. Below are examples.
Standard integration¶
Use esc.integral(F, vp, a, b, numpoints=7, epsabs=None, epsrel=None, maxiter=None) for $\int_a^b F\,d(\text{vp})$, where vp is either a functor variable or a functor parameter. F can be real-valued or complex-valued; a, b are the limits. Options: numpoints ∈ {7, 15, 21, 51, 61}; epsabs and epsrel set absolute and relative tolerances for adaptive refinement; maxiter limits iterations (set to 0 to disable adaptive integration). Return type: if you integrate over the only variable, the result is a dependent parameter; if you integrate over one variable of a multi-variable functor (e.g. $\int F(x,y)\,dx$), the result is a functor of the remaining variable(s). Below are several examples.
X = esc.var("X")
Y = esc.var("Y")
p = esc.par("P", 1.0)
F1 = esc.sin(p * X)
# we create an integral of f(x; p) - I, which is a parameter
I = esc.integral(F1, X, 0, 1.0)
print(type(I))
print(I.num_variables)
# analytical solution: integral_0^1 sin(p*x) dx = (1 - cos(p)) / p
AI = (1 - esc.cos(p)) / p
print(f"Numerical: {I()}")
print(f"Analytical: {AI.value}")
print("Absolute error:", abs(I()-AI.value))
<class 'escape.core.entities.functor_obj'> 0 Numerical: 0.45969769413186035 Analytical: 0.45969769413186023 Absolute error: 1.1102230246251565e-16
# 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)
# analytical solution: integral_0^1 sin(p*X)*cos(p*Y) dX = cos(p*Y) * (1 - cos(p)) / p
AI = (1 - esc.cos(p)) / p * esc.cos(p * Y)
esc.overlay(I, AI).config(
labels=("Numerical", "Analytical"), xlabel="Y", ylabel="I",
title="Integration of 2D functor over X",
)
# integration over parameter: integral_0^1 sin(p*X) dp = (1 - cos(X)) / X
p = esc.par("P", 1.0)
F3 = esc.sin(p * X)
I = esc.integral(F3, p, 0, 1.0, maxiter=100)
AI = (1 - esc.cos(X)) / X
coords = np.linspace(-10, 10, 300)
esc.overlay(I, AI, coordinates=coords).config(
labels=("Numerical", "Analytical"), xlabel="X", ylabel="I",
title="Integration over parameter",
)
Complex integration examples¶
The same esc.integral(...) API supports cplx_functor integrands. Below we use two analytically-known results and compare real and imaginary parts separately with overlays.
- Integration over variable
[ I(Y)=\int_0^{\pi} e^{iX}e^{-Y},dX = 2i,e^{-Y} ]
- Integration over parameter
[ I(X)=\int_0^1 e^{ipX},dp = \frac{e^{iX}-1}{iX} ]
For each case we overlay numerical vs analytical for both real(.) and imag(.).
# --- Complex integration over variable ---
Xc = esc.var("Xc")
Yc = esc.var("Yc")
Fv = esc.exp(1j * Xc) * esc.exp(-Yc)
Iv = esc.integral(Fv, Xc, 0.0, np.pi, numpoints=21, epsabs=1e-10, epsrel=1e-10, maxiter=200)
Av = 2j * esc.exp(-Yc)
coords_y = np.linspace(0.0, 5.0, 300)
esc.overlay(esc.real(Iv), esc.real(Av), esc.imag(Iv), esc.imag(Av), coordinates=coords_y).config(
labels=("Numerical Re", "Analytical Re", "Numerical Im", "Analytical Im"),
xlabel="Y",
ylabel="I",
title="Complex integral over variable: real part",
)
# --- Complex integration over parameter ---
Xp = esc.var("Xp")
Pint = esc.par("Pint", 0.5)
Fp = esc.exp(Pint * Xp* 1j)
Ip = esc.integral(Fp, Pint, 0.0, 1.0, numpoints=21, epsabs=1e-10, epsrel=1e-10, maxiter=200)
Ap = (esc.exp(1j * Xp) - 1.0) / (1j * Xp)
coords_x = np.linspace(0.1, 10.0, 400)
esc.overlay(esc.real(Ip), esc.real(Ap), esc.imag(Ip), esc.imag(Ap), coordinates=coords_x).config(
labels=("Numerical Re", "Analytical Re", "Numerical Im", "Analytical Im"),
xlabel="X",
ylabel="Re(I)",
title="Complex integral over parameter: real part",
)
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.
# Parameter 'p' and its mean value 'p0'
p = esc.par("p", 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$
# coordinates array
x = np.arange(0, 10, 0.01)
G = esc.gamma(X, p, sigma)
WG = G.show(coordinates=x).config(title="Gamma distribution")
# integration over variable
WF = F.show(coordinates=x).config(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, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
WI = I.show(coordinates=x).config(
title="Integration over variable", xlabel="X0", ylabel="I"
)
# integration over parameter
I = esc.average_gamma(F, sigma, p, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10)
WIP = I.show(coordinates=x).config(
title="Integration over parameter", xlabel="X", ylabel="I"
)
esc.show(WG, WF, WI, WIP)
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$
# coordinates array
x = np.arange(0, 10, 0.01)
G = esc.schulz(X, p, sigma)
WG = G.show(coordinates=x).config(title="Schulz-Zimm distribution")
# integration over variable
WF = F.show(coordinates=x).config(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, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
WI = I.show(coordinates=x).config(
title="Integration over variable", xlabel="X0", ylabel="I"
)
# integration over parameter
I = esc.average_schulz(
F, sigma, p, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
WIP = I.show(coordinates=x).config(
title="Integration over parameter", xlabel="X", ylabel="I"
)
esc.show(WG, WF, WI, WIP)
LogNorm distribution function¶
$G(x; x_0, \sigma)=\exp(\log^2(x / x_0) / 2 / \sigma^2))/(\sqrt{2\pi}\sigma x)$
# coordinates array
x = np.arange(0, 10, 0.01)
G = esc.lognorm(X, p, sigma)
WG = G.show(coordinates=x).config(title="LogNorm distribution")
# integration over variable
WF = F.show(coordinates=x).config(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, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
WI = I.show(coordinates=x).config(
title="Integration over variable", xlabel="X0", ylabel="I"
)
# integration over parameter
I = esc.average_lognorm(
F, sigma, p, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
WIP = I.show(coordinates=x).config(
title="Integration over parameter", xlabel="X", ylabel="I"
)
esc.show(WG, WF, WI, WIP)
Normal distribution function¶
$G(x; x_0, \sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{(-(x - x_0)^2 / (2 \sigma^2))}$
# coordinates array
x = np.arange(0, 10, 0.01)
G = esc.normal(X, p, sigma)
WG = G.show(coordinates=x).config(title="Normal distribution")
WF = F.show(coordinates=x).config(title="F(x)=sin(px)^2")
esc.show(WG, WF)
# --- integration over variable ---
I_var = esc.average_normal(
F, sigma * 2.355, X, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
# analytical: <sin^2(px)>_N = 0.5*(1 - cos(2*p*x0)*exp(-2*p^2*sigma^2))
A_var = 0.5 * (1 - esc.cos(2 * p * X) * esc.exp(-2 * p * p * sigma * sigma))
esc.overlay(I_var, A_var, coordinates=x).config(
labels=("Numerical", "Analytical"), xlabel="X0", ylabel="I",
title="Normal: integration over variable",
)
# --- integration over parameter ---
I_par = esc.average_normal(
F, sigma * 2.355, p, maxiter=100, epsrel=1e-8, epsabs=1e-8, numstd=10
)
# analytical: <sin^2(p*x)>_N = 0.5*(1 - cos(2*p0*x)*exp(-2*x^2*sigma^2))
A_par = 0.5 * (1 - esc.cos(2 * p * X) * esc.exp(-2 * X * X * sigma * sigma))
esc.overlay(I_par, A_par, coordinates=x).config(
labels=("Numerical", "Analytical"), xlabel="X", ylabel="I",
title="Normal: integration over parameter",
)
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$
# coordinates array
x = np.arange(0, 10, 0.01)
G = esc.triangular( X, p, sigma)
WG = G.show(coordinates=x).config(title="Triangular distribution")
# integration over variable
WF = F.show(coordinates=x).config(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, maxiter=100, epsrel=1e-8, epsabs=1e-8)
WI = I.show(coordinates=x).config(
title="Integration over variable", xlabel="X0", ylabel="I"
)
# integration over parameter
I = esc.average_triangular(F, sigma * 2, p, maxiter=100, epsrel=1e-8, epsabs=1e-8)
WIP = I.show(coordinates=x).config(
title="Integration over parameter", xlabel="X", ylabel="I"
)
esc.show(WG, WF, WI, WIP)
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$
# coordinates array
x = np.arange(0, 10, 0.01)
w = sigma * 2 # FWHM of the uniform distribution
G = esc.uniform(X, p, sigma)
WG = G.show(coordinates=x).config(title="Uniform distribution")
WF = F.show(coordinates=x).config(title="F(x)=sin(px)^2")
esc.show(WG, WF)
# --- integration over variable ---
I_var = esc.average_uniform(F, w, X, maxiter=100, epsrel=1e-8, epsabs=1e-8)
# analytical: <sin^2(px)>_U = 0.5 - sin(p*w)/(2*p*w) * cos(2*p*x0)
A_var = 0.5 - esc.sin(p * w) / (2 * p * w) * esc.cos(2 * p * X)
esc.overlay(I_var, A_var, coordinates=x).config(
labels=("Numerical", "Analytical"), xlabel="X0", ylabel="I",
title="Uniform: integration over variable",
)
# --- integration over parameter ---
I_par = esc.average_uniform(F, w, p, maxiter=100, epsrel=1e-8, epsabs=1e-8)
# analytical: <sin^2(p*x)>_U = 0.5 - sin(w*x)/(2*w*x) * cos(2*p0*x)
A_par = 0.5 - esc.sin(w * X) / (2 * w * X) * esc.cos(2 * p * X)
esc.overlay(I_par, A_par, coordinates=x).config(
labels=("Numerical", "Analytical"), xlabel="X", ylabel="I",
title="Uniform: integration over parameter",
)
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.
# averaging over variable
X0 = esc.var("X0")
G = 1 / np.sqrt(2 * np.pi) / sigma * esc.exp(-(X-X0 )**2 / (2 * sigma * sigma))
G.variables
[variable(name='X'), variable(name='X0')]
I = esc.average(F, G, X, X0, X - 10 * sigma, X + 10 * sigma, maxiter=100)
# analytical: <sin^2(px)>_gauss = 0.5*(1 - cos(2*p*x0)*exp(-2*p^2*sigma^2))
AI_var = 0.5 * (1 - esc.cos(2 * p * X0) * esc.exp(-2 * p * p * sigma * sigma))
coords = np.arange(0, 10, 0.01)
esc.overlay(I, AI_var, coordinates=coords).config(
labels=("Numerical", "Analytical"), xlabel="X0", ylabel="I",
title="Custom Gaussian: averaging over variable",
)
# averaging over parameter
p0 = esc.par("Mean", 1)
G = 1 / np.sqrt(2 * np.pi) / sigma * esc.exp(-(p - p0) * (p - p0) / (2 * sigma * sigma))
I = esc.average(F, G, p, p0, p - 5 * sigma, p + 5 * sigma, maxiter=100)
# analytical: <sin^2(p*x)>_gauss = 0.5*(1 - cos(2*p0*x)*exp(-2*x^2*sigma^2))
AI_par = 0.5 * (1 - esc.cos(2 * p * X) * esc.exp(-2 * X * X * sigma * sigma))
coords = np.linspace(-10, 10, 10000)
esc.overlay(I, AI_par, coordinates=coords).config(
labels=("Numerical", "Analytical"), xlabel="X", ylabel="I",
title="Custom Gaussian: averaging over parameter",
)