From 880ab4381012aa53f423071da0bc4f10dd9a159a Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 27 Aug 2019 17:36:01 -0500 Subject: [PATCH 1/9] Add an algorithm to compute Clenshaw-Curtis and Fejer's rules --- modepy/quadrature/clenshaw_curtis.py | 102 +++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 modepy/quadrature/clenshaw_curtis.py diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py new file mode 100644 index 0000000..f9dde06 --- /dev/null +++ b/modepy/quadrature/clenshaw_curtis.py @@ -0,0 +1,102 @@ +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): + r"""Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadratures + by DFTs, assuming n >= 2. + + 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 + + """ + N = np.arange(1, n, 2) # noqa + l = len(N) # noqa + m = n - l + K = np.arange(0, m) # noqa + + # Fejer2 nodes: k=0,1,...,n; weights: wf2, wf2_n=wf2_0=0 + v0 = np.concatenate([2 / N / (N - 2), 1 / N[-1:], np.zeros(m)]) + v2 = 0 - v0[:-1] - v0[-1:0:-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) + print('\nFejer-2:') + print(xf2) + print(wf2) + print(np.sum(wf2) - 2) + + # Clenshaw-Curtis nodes: k=0,1,...,n; weights: wcc, wcc_n=wcc_0 + 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) + print('\nClenshaw-Curtis:') + print(xcc) + print(wcc) + print(np.sum(wcc) - 2) + + # 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) + print('\nFejer-1:') + print(xf1) + print(wf1) + print(np.sum(wf1) - 2) + + # tested against ApproxFun's 10th order rules + # https://github.com/JuliaApproximation/ApproxFun.jl/issues/316 + + +class ClenshawCurtisQuadrature(Quadrature): + r"""Clenshaw-Curtis quadrature of order *N*. + + Inherits from :class:`modepy.Quadrature`. See there for the interface + to obtain nodes and weights. + """ + pass + + +if __name__ == "__main__": + _fejer(9) -- GitLab From fc1286a4a8c0b45fe3f727a8f4c25c4248581b91 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 28 Aug 2019 14:36:01 -0500 Subject: [PATCH 2/9] flake8 fixes --- modepy/quadrature/jacobi_gauss.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py index 2235983..2880a65 100644 --- a/modepy/quadrature/jacobi_gauss.py +++ b/modepy/quadrature/jacobi_gauss.py @@ -35,7 +35,7 @@ class JacobiGaussQuadrature(Quadrature): where :math:`\alpha, \beta > -1`. :arg backend: Either ``"builtin"`` or ``"scipy"``. - + When the ``"builtin"`` back-end is in use, there is an additional requirement that :math:`(\alpha, \beta) \not \in \{(-1/2, -1/2)\}`. The ``"scipy"`` backend has no such restriction. @@ -56,7 +56,7 @@ class JacobiGaussQuadrature(Quadrature): backend = 'builtin' if backend == 'builtin': x, w = self.compute_weights_and_nodes(N, alpha, beta) - elif backend== 'scipy': + elif backend == 'scipy': from scipy.special.orthogonal import roots_jacobi x, w = roots_jacobi(N + 1, alpha, beta) else: @@ -186,7 +186,10 @@ def jacobi_gauss_lobatto_nodes(alpha, beta, N, backend=None): # noqa if N == 1: return x - x[1:-1] = np.array(JacobiGaussQuadrature(alpha + 1, beta + 1, N - 2, backend).nodes).real + x[1:-1] = np.array( + JacobiGaussQuadrature( + alpha + 1, beta + 1, N - 2, backend + ).nodes).real return x -- GitLab From 172152687eb8473fc2793357b4a89cbb426efc36 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 28 Aug 2019 15:02:30 -0500 Subject: [PATCH 3/9] Add CC, F1 and F2 quadrature rules --- doc/quadrature.rst | 11 ++ modepy/__init__.py | 2 + modepy/quadrature/clenshaw_curtis.py | 144 ++++++++++++++++++--------- test/test_quadrature.py | 31 ++++++ 4 files changed, 139 insertions(+), 49 deletions(-) diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 4200ff0..46d7d7b 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 02c7298..5e6003f 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 diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py index f9dde06..bd9a023 100644 --- a/modepy/quadrature/clenshaw_curtis.py +++ b/modepy/quadrature/clenshaw_curtis.py @@ -27,9 +27,9 @@ import numpy as np from modepy.quadrature import Quadrature -def _fejer(n): - r"""Weights of the Fejer2, Clenshaw-Curtis and Fejer1 quadratures - by DFTs, assuming n >= 2. +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). @@ -39,64 +39,110 @@ def _fejer(n): 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 - # Fejer2 nodes: k=0,1,...,n; weights: wf2, wf2_n=wf2_0=0 - v0 = np.concatenate([2 / N / (N - 2), 1 / N[-1:], np.zeros(m)]) - v2 = 0 - v0[:-1] - v0[-1:0:-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) - print('\nFejer-2:') - print(xf2) - print(wf2) - print(np.sum(wf2) - 2) - - # Clenshaw-Curtis nodes: k=0,1,...,n; weights: wcc, wcc_n=wcc_0 - 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) - print('\nClenshaw-Curtis:') - print(xcc) - print(wcc) - print(np.sum(wcc) - 2) - - # 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) - print('\nFejer-1:') - print(xf1) - print(wf1) - print(np.sum(wf1) - 2) - - # tested against ApproxFun's 10th order rules - # https://github.com/JuliaApproximation/ApproxFun.jl/issues/316 + 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*. + 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`. """ - pass + 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. -if __name__ == "__main__": - _fejer(9) + 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") + + Quadrature.__init__(self, x, w) diff --git a/test/test_quadrature.py b/test/test_quadrature.py index fd14ce0..865ba6b 100644 --- a/test/test_quadrature.py +++ b/test/test_quadrature.py @@ -71,6 +71,37 @@ 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), -- GitLab From bf7546f58f00c2f3084c6ec8f6f0d9dabec214f9 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 28 Aug 2019 15:04:07 -0500 Subject: [PATCH 4/9] flake8 fix --- test/test_quadrature.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_quadrature.py b/test/test_quadrature.py index 865ba6b..b8f728d 100644 --- a/test/test_quadrature.py +++ b/test/test_quadrature.py @@ -93,6 +93,7 @@ def test_fejer_quadrature(kind): for deg in range(1, 9 + 1): s = deg * 3 quad = FejerQuadrature(s, kind) + def f(x): return x**deg -- GitLab From 39d7a1df93eaa5409154cdf6ea4c1e27e06de06c Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 28 Aug 2019 15:06:49 -0500 Subject: [PATCH 5/9] Add new rules to __all__ --- modepy/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/modepy/__init__.py b/modepy/__init__.py index 5e6003f..1fc1cfb 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -67,6 +67,7 @@ __all__ = [ "GaussLegendreQuadrature", "ChebyshevGaussQuadrature", "XiaoGimbutasSimplexQuadrature", "GrundmannMoellerSimplexQuadrature", "VioreanuRokhlinSimplexQuadrature", + "ClenshawCurtisQuadrature", "FejerQuadrature", ] from pytools import MovedFunctionDeprecationWrapper -- GitLab From b3e033d1245ee9170fddd198b21a99a737cbefa2 Mon Sep 17 00:00:00 2001 From: Xiaoyu Wei Date: Wed, 4 Sep 2019 22:47:13 +0200 Subject: [PATCH 6/9] Apply suggestion to modepy/quadrature/clenshaw_curtis.py --- modepy/quadrature/clenshaw_curtis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py index bd9a023..3c23ac0 100644 --- a/modepy/quadrature/clenshaw_curtis.py +++ b/modepy/quadrature/clenshaw_curtis.py @@ -145,4 +145,4 @@ class FejerQuadrature(Quadrature): else: raise ValueError("kind must be either 1 or 2") - Quadrature.__init__(self, x, w) + super(FejerQuadrature, self).__init__(x, w) -- GitLab From 61ca2e057e809c8e9dd371e2225e9d65e0919e70 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 4 Sep 2019 15:48:08 -0500 Subject: [PATCH 7/9] Fejer raises exception on accessing exact_to --- modepy/quadrature/clenshaw_curtis.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py index bd9a023..33026c7 100644 --- a/modepy/quadrature/clenshaw_curtis.py +++ b/modepy/quadrature/clenshaw_curtis.py @@ -146,3 +146,8 @@ class FejerQuadrature(Quadrature): raise ValueError("kind must be either 1 or 2") Quadrature.__init__(self, x, w) + + @property + def exact_to(self): + raise ValueError("%s has no known exact_to information" + % str(self.__class__.__name__)) -- GitLab From 48b704b4108534fda5e1cc60064dc1392908de17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 4 Sep 2019 23:14:18 +0200 Subject: [PATCH 8/9] Apply suggestion to modepy/quadrature/clenshaw_curtis.py --- modepy/quadrature/clenshaw_curtis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py index 41a29cb..f8521de 100644 --- a/modepy/quadrature/clenshaw_curtis.py +++ b/modepy/quadrature/clenshaw_curtis.py @@ -150,4 +150,4 @@ class FejerQuadrature(Quadrature): @property def exact_to(self): raise ValueError("%s has no known exact_to information" - % str(self.__class__.__name__)) + % type(self).__name__)) -- GitLab From b1d9c37d06997ec32fa90ab8e5fbf7f7550cda6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Wed, 4 Sep 2019 23:15:42 +0200 Subject: [PATCH 9/9] Apply suggestion to modepy/quadrature/clenshaw_curtis.py --- modepy/quadrature/clenshaw_curtis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py index f8521de..f8d5b7f 100644 --- a/modepy/quadrature/clenshaw_curtis.py +++ b/modepy/quadrature/clenshaw_curtis.py @@ -150,4 +150,4 @@ class FejerQuadrature(Quadrature): @property def exact_to(self): raise ValueError("%s has no known exact_to information" - % type(self).__name__)) + % type(self).__name__) -- GitLab