psf_generator.utils.integrate#

A collection of numerical integration rules in 1D.

We provide two common rules written in PyTorch: trapezoid and simpson, equivalent to their counterparts in scipy.integrate.

How it works is briefly explained here:

We would like to compute the integral \(I = \int_{a}^{b} f(x) dx\) of a function \(f(x)\) over an interval \([a, b]\). We partition \([a, b]\) into \(n\) equidistant sub-intervals by \(N=n+1\) nots \(a \leq x_0 < x_1 < \ldots < x_n \leq b\), in other words, \(x_k = a + kh, k=0, \ldots, n\) with stepsize \(h=\frac{b - a}{n}\).

The composite trapezoid rule

\[T_n = \frac{h}{2}\left(f(a) + 2\sum_{k=1}^{n-1}f(x_k) + f(b)\right).\]

The composite Simpson’s rule

\[S_n = \frac{h}{6}\left(f(a) + 4\sum_{k=1}^{n-1}f{x_{k+\frac{1}{2}}} + 2\sum_{k=1}^{n-1}f(x_k) + f(b)\right),\]

where \(x_{k+\frac{1}{2}} = \frac{x_k + x_{k+1}}{2}\).

In implementation, trapezoid is written exactly as its formula. Simpson’s rule requires the function value of the midpoint which is not provided, we view the partition stepsize as \(h = 2\frac{b - a}{n}\) and the odd nots as the midpoint instead. The formula then becomes

\[S_n = \frac{h}{3}\left(f(a) + 4\sum_{k=1}^{n/2}f(x_{2k-1}) + 2\sum_{k=1}^{n/2-1}f(x_{2k}) + f(b)\right).\]

Functions#

riemann_rule(→ torch.Tensor)

Riemann quadrature rule of precision \(O(h)\).

simpsons_rule(→ torch.Tensor)

Composite Simpson's rule, see also [2].

Module Contents#

psf_generator.utils.integrate.riemann_rule(fs: torch.Tensor, dx: float) → torch.Tensor[source]#

Riemann quadrature rule of precision \(O(h)\).

Parameters#

fstorch.Tensor

The integrand evaluations of shape (N, number_of_integrals).

dxfloat

Bin width or step size for evaluation \(h = 1 / (N - 1)\).

Returns#

output: torch.Tensor

Integral evaluated by Riemann rule of shape (num_integrals,).

psf_generator.utils.integrate.simpsons_rule(fs: torch.Tensor, dx: float) → torch.Tensor[source]#

Composite Simpson’s rule, see also [2].

Parameters#

fstorch.Tensor

The integrand evaluations of shape (N, number_of_integrals).

dxfloat

Step size.

Returns#

output: torch.Tensor

Integral evaluated by Simpson’s rule of shape (num_integrals,).

Notes#

  • \(h = \frac{b - a}{N - 1}\).

  • Simpson’s rule only works correctly with grids of odd sizes (i.e. \(N = 2^K + 1\)).

References#

On this page