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).