# -*- coding: utf-8 -*-
"""Let :math:`f(\\xi)` be an integrable function defined in :math:`[-1,1]`.
The integral of :math:`f(\\xi)` over :math:`[-1,1]` can be
approximated by
.. math::
\int_{-1}^{1}f(\\xi)\\approx\\sum_{i} f(\\xi_{i})\\omega_{i},
where :math:`\\xi_{i}` and :math:`\\omega_{i}` are the quadrature nodes
and weights. For example, degree 3 Gauss quadrature nodes are around
.. math::
[-0.86,\ -0.34,\ 0.34,\ 0.86],
and degree 3 Gauss quadrature weights are around
.. math::
[0.35,\ 0.65,\ 0.65,\ 0.35].
The integral of function :math:`f(\\xi) = \\sin(\\pi\\xi)` over
:math:`[-1,1]` can be approximated by
.. math::
\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 :math:`f(\\xi,\\eta,\\varsigma)` be an
integrable function in :math:`\\Omega_{\\mathrm{ref}} = [-1,1]^3`.
Its integral over :math:`\\Omega_{\\mathrm{ref}}` can be computed
numerically through
.. math::
\\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 :math:`\\xi_{i}` and :math:`\\omega_{i}`,
:math:`\\eta_{j}` and :math:`\\omega_{j}` and
:math:`\\varsigma_{k}` and :math:`\\omega_{k}` are three groups of
1D quadrature nodes and weights.
For an integral of :math:`f(x,y,z)` over a domain other
than :math:`\\Omega_{\\mathrm{ref}}`, for example, over
:math:`\\Omega = \\Phi(\\Omega_{\\mathrm{ref}})`, the numerical integral
is
.. math::
\\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 :math:`(x_i, y_j, z_k) = \Phi(\\xi_{i},\\eta_{j},\\varsigma_{k})` and
:math:`g` is the metric of the mapping :math:`\\Phi`, see :ref:`ct`.
For example, if we integrate
:math:`f(x,y,z) = \\pi^3\\sin(\\pi x)\\sin(\\pi y)\\sin(\\pi z)` over domain
:math:`[0, 1]^3`, we know analytically that
:math:`\\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
:download:`[quadrature.py]</contents/LIBRARY/ptc/mathischeap_ptc/quadrature.py>`.
Dependence may exist. In case of error, check import and install required
packages or download required scripts. © mathischeap.com
"""
import numpy as np
from scipy.special import legendre
def _newton_method_(f, df_dx, x_0, n_max):
"""Newton method for root finding."""
x = [x_0,]
for i in range(n_max - 1):
# noinspection PyUnresolvedReferences
x.append(x[i] - f(x[i]) / df_dx(x[i]))
if abs(x[i + 1] - x[i]) < np.finfo(float).eps * 10:
break
return x[-1]
[docs]
def Lobatto(p):
"""Compute the Lobatto quadrature nodes and weights of degree ``p``
using newton method. Credits belong to Lorenzo.
:param p: The quadrature degree (positive integer).
:type p: int
:return: 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])
"""
x_0 = np.cos(np.arange(1, p) / p * np.pi)
nodal_pts = np.zeros((p + 1))
nodal_pts[0] = 1
nodal_pts[-1] = -1
legendre_prime = lambda x, n: (n * legendre(n - 1)(x) -
n * x * legendre(n)(x)) / (1 - x ** 2)
for i, ch_pt in enumerate(x_0):
leg_p = lambda x: (1 - x**2)**2 * legendre_prime(x,p)
leg_pp = lambda x: (2 * x * legendre_prime(x, p) -
p * (p + 1) * legendre(p)(x)) * (1 - x ** 2)
nodal_pts[i + 1] = _newton_method_(leg_p, leg_pp, ch_pt, 100)
weights = 2 / (p * (p + 1) * (legendre(p)(nodal_pts))**2)
return nodal_pts[::-1], weights
[docs]
def Gauss(p):
"""Compute the Gauss quadrature nodes and weights of degree ``p``
with the numpy package.
:param p: The quadrature degree (positive integer).
:type p: int
:return: 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])
"""
return np.polynomial.legendre.leggauss(p+1)
if __name__ == "__main__":
import doctest
doctest.testmod()
