From 0d9bf863a9269b318c3538126e3c9b2d34007b09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Wed, 30 Sep 2026 14:59:27 -0400 Subject: [PATCH 1/3] Fix FejerFirst and FejerSecond implementations --- src/grid/onedgrid.py | 47 ++++++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/src/grid/onedgrid.py b/src/grid/onedgrid.py index 0ca8eeb8..939261f3 100644 --- a/src/grid/onedgrid.py +++ b/src/grid/onedgrid.py @@ -133,7 +133,7 @@ class GaussChebyshev(OneDGrid): x_i =& \cos\left( \frac{2i-1}{2n}\pi \right) \\ w_i =& \frac{\pi}{n} - However, to integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: + To integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: .. math:: \int_{-1}^{1}g(x)dx \approx \sum_{i=1}^n \left(w_i\sqrt{1-x_i^2}\right)g(x_i) = @@ -212,7 +212,7 @@ class GaussChebyshevType2(OneDGrid): x_i =& \cos\left( \frac{i}{n+1} \pi \right) \\ w_i =& \frac{\pi}{n+1} \sin^2 \left( \frac{i}{n+1} \pi \right) - However, to integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: + To integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: .. math:: \int_{-1}^{1} g(x) dx \approx \sum_{i=1}^n \left(\frac{w_i}{\sqrt{1-x_i^2}}\right) g(x_i) = @@ -256,7 +256,7 @@ class GaussChebyshevLobatto(OneDGrid): w_{1} = w_{n} =& \frac{\pi}{2(n-1)} \\ w_{i\neq 1,n} =& \frac{\pi}{n-1} - However, to integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: + To integrate a given function :math:`g(x)` over :math:`[-1, 1]`\, this is re-written as: .. math:: \int_{-1}^{1}g(x) dx \approx \sum_{i=1}^n \left(w_i \sqrt{1-x_i^2}\right) g(x_i) = @@ -671,9 +671,9 @@ class FejerFirst(OneDGrid): w_i &= \frac{2}{n}\bigg(1 - 2 \sum_{j=1}^{\lfloor n/2 \rfloor} \frac{\cos(2j \theta_j)}{4 j^2 - 1} \bigg), - where :math:`k=1,\cdots, n`\. It uses the zeros of the Chebyshev polynomial. - If discontinuous, it is recommended to break the intervals at the discontinuities - and handled separately. + where :math:`i=1,\cdots, n` and :math:`n` is the number of quadrature points. It uses the zeros + of the Chebyshev polynomial. If discontinuous, it is recommended to break the intervals at the + discontinuities and handled separately. """ @@ -700,15 +700,16 @@ def __init__(self, npoints: int): points = np.cos(theta) nsum = npoints // 2 - j = np.arange(nsum - 1) + 1 + j = np.arange(nsum) + 1 - bj = 2.0 * np.ones(nsum - 1) / (4 * j**2 - 1) + bj = 2.0 / (4 * j**2 - 1) cij = np.cos(np.outer(2 * j, theta)) di = bj @ cij - weights = 1 - di + weights = (1 - di) * (2 / npoints) + # reverse the order of points and weights to have them in ascending order points = points[::-1] - weights = weights[::-1] * (2 / npoints) + weights = weights[::-1] super().__init__(points, weights, (-1, 1)) @@ -720,15 +721,14 @@ class FejerSecond(OneDGrid): The definition of this quadrature is: .. math:: - \theta_i &= k \pi / n \\ + \theta_i &= \frac{i\pi}{n + 1} \\ x_i &= \cos(\theta_i) \\ - w_i &= \frac{4 \sin(\theta_i)}{n} \sum_{j=1}^{\lfloor n/2 \rfloor} + w_i &= \frac{4 \sin(\theta_i)}{n + 1} \sum_{j=1}^{\lfloor (n + 1)/2 \rfloor} \frac{\sin(2j - 1)\theta_i}{2j - 1}\\ - where :math:`k=1, \cdots n - 1` and :math:`n` is the number of points. This - method is considered more practical than the first method. If discontinuous, it is + where :math:`i = 1, \ldots, n` and :math:`n` is the number of quadrature + points. This method is considered more practical than the first method. If discontinuous, it is recommended to break the intervals at the discontinuities and handled separately. - """ name = "Fejer-Second" @@ -747,22 +747,23 @@ def __init__(self, npoints: int): One-dimensional grid instance. """ if npoints <= 1: - raise ValueError("npoints must be greater that one, given {npoints}") + raise ValueError(f"npoints must be greater than one, given {npoints}") theta = np.pi * (np.arange(npoints) + 1) / (npoints + 1) - points = np.cos(theta) nsum = (npoints + 1) // 2 - j = np.arange(nsum - 1) + 1 + j = np.arange(nsum) + 1 - bj = np.ones(nsum - 1) / (2 * j - 1) - sij = np.sin(np.outer(2 * j - 1, theta)) - wi = bj @ sij - weights = 4 * np.sin(theta) * wi + odd_indices = 2 * j - 1 + sin_coeffs = 1.0 / odd_indices + sin_terms_ij = np.sin(np.outer(odd_indices, theta)) + weighted_sum = sin_coeffs @ sin_terms_ij + weights = 4 * np.sin(theta) / (npoints + 1) * weighted_sum + # sort the points and weights in ascending order points = points[::-1] - weights = weights[::-1] / (npoints + 1) + weights = weights[::-1] super().__init__(points, weights, (-1, 1)) From 712361b59f4bc3743d9f677893be2da82faa16d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Wed, 30 Sep 2026 15:12:45 -0400 Subject: [PATCH 2/3] Fix FejerSecond test The test was a clone of the erroneous implementation. --- src/grid/tests/test_onedgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/tests/test_onedgrid.py b/src/grid/tests/test_onedgrid.py index 641fd619..99645577 100644 --- a/src/grid/tests/test_onedgrid.py +++ b/src/grid/tests/test_onedgrid.py @@ -463,7 +463,7 @@ def test_FejerSecond(self): for k in range(0, 10): serie = 0 - for m in range(1, nsum): + for m in range(1, nsum + 1): serie += np.sin((2 * m - 1) * theta[k]) / (2 * m - 1) weights[k] = (4 * np.sin(theta[k]) / 11) * serie From fb8f2cbaae3e31f848c7322e242c073bf80b53ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Thu, 1 Oct 2026 16:04:38 -0400 Subject: [PATCH 3/3] Add polynomial exactness test and docstring fixes for Fejer1 and Fejer2 --- src/grid/onedgrid.py | 4 ++-- src/grid/tests/test_onedgrid.py | 17 +++++++++++++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/grid/onedgrid.py b/src/grid/onedgrid.py index 939261f3..55acd195 100644 --- a/src/grid/onedgrid.py +++ b/src/grid/onedgrid.py @@ -669,7 +669,7 @@ class FejerFirst(OneDGrid): \theta_i &= \frac{(2i - 1)\pi}{2n}, \\ x_i &= \cos(\theta_i), \\ w_i &= \frac{2}{n}\bigg(1 - 2 \sum_{j=1}^{\lfloor n/2 \rfloor} - \frac{\cos(2j \theta_j)}{4 j^2 - 1} \bigg), + \frac{\cos(2j \theta_i)}{4 j^2 - 1} \bigg), where :math:`i=1,\cdots, n` and :math:`n` is the number of quadrature points. It uses the zeros of the Chebyshev polynomial. If discontinuous, it is recommended to break the intervals at the @@ -724,7 +724,7 @@ class FejerSecond(OneDGrid): \theta_i &= \frac{i\pi}{n + 1} \\ x_i &= \cos(\theta_i) \\ w_i &= \frac{4 \sin(\theta_i)}{n + 1} \sum_{j=1}^{\lfloor (n + 1)/2 \rfloor} - \frac{\sin(2j - 1)\theta_i}{2j - 1}\\ + \frac{\sin((2j - 1)\theta_i)}{2j - 1}\\ where :math:`i = 1, \ldots, n` and :math:`n` is the number of quadrature points. This method is considered more practical than the first method. If discontinuous, it is diff --git a/src/grid/tests/test_onedgrid.py b/src/grid/tests/test_onedgrid.py index 99645577..436cb3f8 100644 --- a/src/grid/tests/test_onedgrid.py +++ b/src/grid/tests/test_onedgrid.py @@ -474,6 +474,23 @@ def test_FejerSecond(self): assert_allclose(grid.points, points) assert_allclose(grid.weights, weights) + def test_fejer_polynomial_exactness(self): + """Test polynomial exactness of Fejer quadrature rules.""" + for grid_class in (FejerFirst, FejerSecond): + for npoints in range(2, 16): + grid = grid_class(npoints) + + for degree in range(npoints): + result = grid.integrate(grid.points**degree) + # even polynomials: integral of x^n over [-1, 1] is 2/(n+1) + if degree % 2 == 0: + expected = 2.0 / (degree + 1) + # odd polynomials: integral of x^n over [-1, 1] is 0 (asymmetric function) + else: + expected = 0.0 + + assert_allclose(result, expected, rtol=1e-12, atol=1e-12) + def test_AuxiliarTrefethenSausage(self): """Test for Auxiliary functions using in Trefethen Sausage.""" xref = np.array([-1.0, 0.0, 1.0])