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
>>> results = list()
...     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

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


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