diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 4200ff0d289d369c715700f7fc7e3aa323827b68..46d7d7bbad81769f0743001b67554c2fc1d052e0 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -35,6 +35,17 @@ Jacobi-Gauss quadrature in one dimension .. autofunction:: legendre_gauss_lobatto_nodes +Clenshaw-Curtis and Fejér quadrature in one dimension +----------------------------------------------------- + +.. automodule:: modepy.quadrature.clenshaw_curtis + +.. currentmodule:: modepy + +.. autoclass:: ClenshawCurtisQuadrature + +.. autoclass:: FejerQuadrature + Quadratures on the simplex -------------------------- diff --git a/modepy/__init__.py b/modepy/__init__.py index 9fa1d3143c6557688c0ac9f0d634d54aaa4da3f4..5931c6f9bcf052b8334a1b27d978c755668e19c6 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -40,6 +40,8 @@ from modepy.quadrature.jacobi_gauss import ( from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature +from modepy.quadrature.clenshaw_curtis import ( + ClenshawCurtisQuadrature, FejerQuadrature) from modepy.version import VERSION_TEXT as __version__ # noqa: N811 @@ -66,6 +68,7 @@ __all__ = [ "GaussGegenbauerQuadrature", "XiaoGimbutasSimplexQuadrature", "GrundmannMoellerSimplexQuadrature", "VioreanuRokhlinSimplexQuadrature", + "ClenshawCurtisQuadrature", "FejerQuadrature", ] from pytools import MovedFunctionDeprecationWrapper diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py new file mode 100644 index 0000000000000000000000000000000000000000..f8d5b7f8db0c819701e5ca9819d0328edbf24cd8 --- /dev/null +++ b/modepy/quadrature/clenshaw_curtis.py @@ -0,0 +1,153 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import numpy as np + +from modepy.quadrature import Quadrature + + +def _fejer(n, rule): + r"""Nodes and weights of the Fejer2, Clenshaw-Curtis and Fejer1 + quadratures by DFTs. + + Nodes: x_k = cos(k * pi / n). + + The algorithm follows: + Jörg Waldvogel, + Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules, + BIT Numerical Mathematics, 2003, Vol. 43, No. 1, pp. 001-018, + https://doi.org/10.1007/s10543-006-0045-4 + + :arg n: the scheme's order. n is related to the number of quadrature + nodes as follows: + - for f1, n_1_nodes = n, n >= 1 + - for f2, n_q_nodes = n - 1, n >= 2 + - for cc, n_q_nodes = n + 1, n >= 1 + (a minimal cc rule must include the end points) + + :arg rule: any of ``"cc"``, ``"f1"`` or ``"f2"``. + """ + N = np.arange(1, n, 2) # noqa + l = len(N) # noqa + m = n - l + K = np.arange(0, m) # noqa + + if not n >= 1: + raise RuntimeError('Invalid n = %d (must be n >= 1)' % n) + + if rule not in ['f1', 'f2', 'cc']: + raise NotImplementedError('Invalide rule: %s' % rule) + + if rule == 'f2' or rule == 'cc': + + v0 = np.concatenate([2 / N / (N - 2), 1 / N[-1:], np.zeros(m)]) + v2 = 0 - v0[:-1] - v0[-1:0:-1] + + if rule == 'f2': + # Fejer2 nodes: k=0,1,...,n; weights: wf2, wf2_n=wf2_0=0 + # nodes with zero weights are not included in the return values + if n == 1: + raise RuntimeError("Fejer1 does not exist for n=1") + wf2 = np.fft.ifft(v2) + assert np.allclose(wf2.imag, 0) + wf2 = wf2.real[1:] + xf2 = np.cos(np.arange(1, n) * np.pi / n) + return xf2, wf2 + + if rule == 'cc': + # Clenshaw-Curtis nodes: k=0,1,...,n; weights: wcc, wcc_n=wcc_0 + if n == 1: + return np.array([-1, 1]), np.array([1, 1]) + g0 = -np.ones(n) + g0[l] = g0[l] + n + g0[m] = g0[m] + n + g = g0 / (n**2 - 1 + (n % 2)) + wcc = np.fft.ifft(v2 + g) + assert np.allclose(wcc.imag, 0) + wcc = np.concatenate([wcc, wcc[:1]]).real + xcc = np.cos(np.arange(0, n + 1) * np.pi / n) + return xcc, wcc + + if rule == 'f1': + # Fejer1 nodes: k=1/2,3/2,...,n-1/2; vector of weights: wf1 + v0 = np.concatenate( + [2 * np.exp(1j * np.pi * K / n) / (1 - 4 * (K**2)), + np.zeros(l+1)]) + v1 = v0[:-1] + np.conj(v0[-1:0:-1]) + wf1 = np.fft.ifft(v1) + assert np.allclose(wf1.imag, 0) + wf1 = wf1.real + xf1 = np.cos((np.arange(0, n) + 0.5) * np.pi / n) + return xf1, wf1 + + +class ClenshawCurtisQuadrature(Quadrature): + r"""Clenshaw-Curtis quadrature of order *N* (having *N + 1* points). + + Inherits from :class:`modepy.Quadrature`. See there for the interface + to obtain nodes and weights. + + Integrates on the interval :math:`(-1, 1)`. + The quadrature rule is exact up to degree :math:`N`; however, its + performance for differentiable functions is comparable with the classic + Gauss-Legendre quadrature which is exact for polynomials of degree up + to :math:`2N + 1`. + """ + def __init__(self, N): # noqa + x, w = _fejer(N, 'cc') + self.exact_to = N + Quadrature.__init__(self, x, w) + + +class FejerQuadrature(Quadrature): + r"""Fejér's quadrature rules of order *N*, categorized in two kinds. + The Fejér's quadrature rule of first kind has *N* points; while the + second kind has *N - 1* points. + + The first kind uses Chebyshev nodes, i.e. roots of the Chebyshev + polynomials. The second kind uses the interior extrema of the Chebyshev + polynomials, i.e. the true stationary points. + + The second-kind Fejér's quadrature rule is nearly identical to + Clenshaw-Curtis. Both can also be nested. + + Inherits from :class:`modepy.Quadrature`. See there for the interface + to obtain nodes and weights. + + Integrates on the interval :math:`(-1, 1)`. + """ + def __init__(self, N, kind=1): # noqa + if kind == 1: + x, w = _fejer(N, 'f1') + elif kind == 2: + x, w = _fejer(N, 'f2') + else: + raise ValueError("kind must be either 1 or 2") + + super(FejerQuadrature, self).__init__(x, w) + + @property + def exact_to(self): + raise ValueError("%s has no known exact_to information" + % type(self).__name__) diff --git a/test/test_quadrature.py b/test/test_quadrature.py index fd14ce0e124c92fc3321a4f1f1de7b8adcdb20ab..b8f728d0a9db9f50a08464d248345f69b415702b 100644 --- a/test/test_quadrature.py +++ b/test/test_quadrature.py @@ -71,6 +71,38 @@ def test_gauss_quadrature(backend): assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) +def test_clenshaw_curtis_quadrature(): + from modepy.quadrature.clenshaw_curtis import ClenshawCurtisQuadrature + + for s in range(1, 9 + 1): + quad = ClenshawCurtisQuadrature(s) + for deg in range(quad.exact_to + 1): + def f(x): + return x**deg + + i_f = quad(f) + i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + err = abs(i_f - i_f_true) + assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) + + +@pytest.mark.parametrize("kind", [1, 2]) +def test_fejer_quadrature(kind): + from modepy.quadrature.clenshaw_curtis import FejerQuadrature + + for deg in range(1, 9 + 1): + s = deg * 3 + quad = FejerQuadrature(s, kind) + + def f(x): + return x**deg + + i_f = quad(f) + i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + err = abs(i_f - i_f_true) + assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) + + @pytest.mark.parametrize(("quad_class", "highest_order"), [ (mp.XiaoGimbutasSimplexQuadrature, None), (mp.VioreanuRokhlinSimplexQuadrature, None),