quadrature.py

Let f(\xi) be an integrable function defined in [-1,1]. The integral of f(\xi) over [-1,1] can be approximated by

\int_{-1}^{1}f(\xi)\approx\sum_{i} f(\xi_{i})\omega_{i},

where \xi_{i} and \omega_{i} are the quadrature nodes and weights. For example, degree 3 Gauss quadrature nodes are around

[-0.86,\ -0.34,\ 0.34,\ 0.86],

and degree 3 Gauss quadrature weights are around

[0.35,\ 0.65,\ 0.65,\ 0.35].

The integral of function f(\xi) = \sin(\pi\xi) over [-1,1] can be approximated by

\int_{-1}^{1}\sin(\pi\xi) \approx 0.35\sin(- 0.86\pi) + 0.65\sin(-0.34\pi) + 0.65\sin(0.34\pi) + 0.35\sin(0.86\pi).

Higher dimensional numerical integral can be computed using the tensor product. For example, in 3D, let f(\xi,\eta,\varsigma) be an integrable function in \Omega_{\mathrm{ref}} = [-1,1]^3. Its integral over \Omega_{\mathrm{ref}} can be computed numerically through

\iiint_{\Omega_{\mathrm{ref}}} f(\xi,\eta,\varsigma) \approx \sum_{i}\sum_{j} \sum_{k} f(\xi_{i},\eta_{j},\varsigma_{k})\omega_{i} \omega_{j}\omega_{k}

where \xi_{i} and \omega_{i}, \eta_{j} and \omega_{j} and \varsigma_{k} and \omega_{k} are three groups of 1D quadrature nodes and weights.

For an integral of f(x,y,z) over a domain other than \Omega_{\mathrm{ref}}, for example, over \Omega = \Phi(\Omega_{\mathrm{ref}}), the numerical integral is

\iiint_{\Omega} f(x,y,z) \approx \sum_{i}\sum_{j} \sum_{k} \sqrt{g(\xi_{i},\eta_{j},\varsigma_{k})} f(x_i, y_j, z_k)\omega_{i} \omega_{j}\omega_{k}

where (x_i, y_j, z_k) = \Phi(\xi_{i},\eta_{j},\varsigma_{k}) and g is the metric of the mapping \Phi, see coordinate_transformation.py.

For example, if we integrate f(x,y,z) = \pi^3\sin(\pi x)\sin(\pi y)\sin(\pi z) over domain [0, 1]^3, we know analytically that \iiint_{[0,1]^3} f(x,y,z) = 8. The following codes do it numerically using Gauss quadrature at degrees 1,2,3,4.

>>> from numpy import pi, sin, meshgrid, tensordot, sum
>>> def f(x,y,z): return sin(pi * x) * sin(pi * y) * sin(pi * z) * pi**3
>>> quad_degrees = [1,2,3,4]
>>> results = list()
>>> for qd in quad_degrees:
...     nodes, weights = Gauss(qd) # use some quadrature degree along 3 directions.
...     mapped_nodes = (nodes + 1) * 0.5
...     xyz = meshgrid(mapped_nodes, mapped_nodes, mapped_nodes, indexing='ij')
...     weights_2d = tensordot(weights, weights, axes = 0)
...     weights_3d = tensordot(weights_2d, weights, axes = 0)
...     numerical_integration = sum(0.125 * f(*xyz) * weights_3d) # g = 0.125 is constant.
...     results.append(numerical_integration)
>>> results
[7.254285290478619, 8.016678540458308, 7.999810742985112, 8.000001323413734]

Obviously, the numerical integration converges to the true value, 8.

Note that we do not wrap the numerical integration into a general function because we will (and also hope you could) use it freely in different circumstances. Here we only provide the functions to compute the quadrature nodes and weights. Two types of them, Gauss quadrature and Lobatto quadrature, are provided.

⭕ To access the source code, click on the [source] button at the right side or click on [quadrature.py]. Dependence may exist. In case of error, check import and install required packages or download required scripts. © mathischeap.com

quadrature.Gauss(p)[source]

Compute the Gauss quadrature nodes and weights of degree p with the numpy package.

Parameters:

p (int) – The quadrature degree (positive integer).

Returns:

A tuple of two outputs:

  1. (numpy.ndarray) - The Gauss quadrature nodes.

  2. (numpy.ndarray) - The Gauss quadrature weights.

Example:

>>> n, w = Gauss(3)
>>> n
array([-0.86113631, -0.33998104,  0.33998104,  0.86113631])
>>> w
array([0.34785485, 0.65214515, 0.65214515, 0.34785485])
quadrature.Lobatto(p)[source]

Compute the Lobatto quadrature nodes and weights of degree p using newton method. Credits belong to Lorenzo.

Parameters:

p (int) – The quadrature degree (positive integer).

Returns:

A tuple of two outputs:

  1. (numpy.ndarray) - The Lobatto quadrature nodes.

  2. (numpy.ndarray) - The Lobatto quadrature weights.

Example:

>>> n, w = Lobatto(3)
>>> n
array([-1.       , -0.4472136,  0.4472136,  1.       ])
>>> w
array([0.16666667, 0.83333333, 0.83333333, 0.16666667])

↩️ Back to Ph.D. thesis complements (ptc).