quadrature.py¶
Let be an integrable function defined in . The integral of over can be approximated by
where and are the quadrature nodes and weights. For example, degree 3 Gauss quadrature nodes are around
- and degree 3 Gauss quadrature weights are around
The integral of function over can be approximated by
Higher dimensional numerical integral can be computed using the tensor product. For example, in 3D, let be an integrable function in . Its integral over can be computed numerically through
where and , and and and are three groups of 1D quadrature nodes and weights.
For an integral of over a domain other than , for example, over , the numerical integral is
where and is the metric of the mapping , see coordinate_transformation.py.
For example, if we integrate over domain , we know analytically that . 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:
(numpy.ndarray) - The Gauss quadrature nodes.
(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:
(numpy.ndarray) - The Lobatto quadrature nodes.
(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).