Fix FejerFirst and FejerSecond implementations

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) =
@@ -669,11 +669,11 @@ 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:`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}
-            \frac{\sin(2j - 1)\theta_i}{2j - 1}\\
+        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))
 

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