Fix accuracy-order loss for even derivatives on non-uniform grids

CHANGES.md

@@ -7,6 +7,7 @@
 - Iterative solvers for PDEs (CG, GMRES, BiCGSTAB, LGMRES, MINRES) with preconditioner support
 - Eigenvalue problem support (eigs and eigsh methods)
 - Robin (mixed) boundary conditions for PDEs
+- Fix order-of-accuracy loss for even-order derivatives on non-uniform grids (central stencil was undersized)
 
 ## Version 0.13.x
 

findiff/coefs.py

@@ -159,13 +159,16 @@ def coefficients_non_uni(deriv, acc, coords, idx):
     _validate_acc(acc)
 
     num_central = 2 * math.floor((deriv + 1) / 2) - 1 + acc
-    num_side = num_central // 2
 
     if deriv % 2 == 0:
         num_coef = num_central + 1
     else:
         num_coef = num_central
 
+    # Non-uniform central stencils get no symmetry bonus, so size them like the
+    # one-sided schemes (num_coef points) to reach the requested order `acc`.
+    num_side = num_coef // 2
+
     if idx < num_side:
         matrix = _build_matrix_non_uniform(0, num_coef - 1, coords, idx)
 
@@ -219,8 +222,10 @@ def calc_coefs_non_uni_batched(deriv, acc, coords):
     _validate_acc(acc)
 
     num_central = 2 * math.floor((deriv + 1) / 2) - 1 + acc
-    num_side = num_central // 2
     num_coef = num_central + 1 if deriv % 2 == 0 else num_central
+    # Non-uniform central stencils get no symmetry bonus, so size them like the
+    # one-sided schemes (num_coef points) to reach the requested order `acc`.
+    num_side = num_coef // 2
 
     N = len(coords)
 

tests/test_coefs.py

@@ -170,12 +170,18 @@ def test_non_uniform(analytic_inv):
     x = np.linspace(0, 10, 100)
     dx = x[1] - x[0]
 
-    c_uni = coefficients(deriv=2, acc=2, analytic_inv=analytic_inv)
-    coefs_uni = c_uni["center"]["coefficients"] / dx**2
-
     c_non_uni = coefficients_non_uni(deriv=2, acc=2, coords=x, idx=5)
     coefs_non_uni = c_non_uni["coefficients"]
 
+    # On a uniform grid the non-uniform scheme must reproduce the equidistant
+    # coefficients for the stencil it actually uses.
+    coefs_uni = (
+        calc_coefs(2, list(c_non_uni["offsets"]), analytic_inv=analytic_inv)[
+            "coefficients"
+        ]
+        / dx**2
+    )
+
     np.testing.assert_array_almost_equal(coefs_non_uni, coefs_uni)
 
 

tests/test_nonuniform_accuracy.py

@@ -0,0 +1,86 @@
+"""Accuracy-order regression tests for non-uniform grids.
+
+On a non-uniform grid the central stencil earns no symmetry bonus, so the
+interior scheme for an even-order derivative must carry as many points as the
+one-sided schemes to reach the requested order. Before that fix the even-order
+interior schemes were one point short and only reached order ``acc - 1`` (the
+non-uniform Laplacian was only first-order in the interior). Odd orders were
+already correct. These tests pin the order down by polynomial exactness and by
+a graded-grid convergence study.
+"""
+
+import numpy as np
+import pytest
+
+from findiff import Diff
+
+
+def _graded_grid(n, ratio=2.0, length=3.0):
+    """Grid with alternating spacings h, ratio*h, ... (fixed O(1) asymmetry).
+
+    A smoothly mapped grid hides the order loss (its asymmetry vanishes as the
+    grid is refined); a graded/boundary-layer mesh does not.
+    """
+    steps = np.empty(n - 1)
+    steps[0::2] = 1.0
+    steps[1::2] = ratio
+    x = np.concatenate([[0.0], np.cumsum(steps)])
+    return x / x[-1] * length
+
+
+@pytest.mark.parametrize("acc", [2, 4])
+@pytest.mark.parametrize("deriv", [1, 2, 3, 4])
+def test_nonuniform_polynomial_exactness(deriv, acc):
+    """An order-acc scheme is exact for polynomials of degree deriv+acc-1 at
+    every grid point. This failed for even derivatives before the fix."""
+    x = _graded_grid(40)
+    degree = deriv + acc - 1
+    p = np.poly1d(np.random.default_rng(0).standard_normal(degree + 1))
+    f = p(x)
+    op = Diff(0, x, acc=acc) ** deriv if deriv > 1 else Diff(0, x, acc=acc)
+    got = op(f)
+    exact = p.deriv(deriv)(x)
+    scale = max(1.0, np.max(np.abs(exact)))
+    assert np.max(np.abs(got - exact)) / scale < 1e-8
+
+
+@pytest.mark.parametrize("deriv,acc", [(1, 2), (1, 4), (2, 2), (2, 4)])
+def test_nonuniform_convergence_order(deriv, acc):
+    """Interior convergence order on a graded grid reaches the requested acc.
+
+    Even derivatives dropped to order acc-1 before the central stencil was
+    widened; odd derivatives were already correct.
+    """
+
+    def f_and_deriv(x):
+        f = np.exp(np.sin(2 * x))
+        if deriv == 1:
+            return f, 2 * np.cos(2 * x) * f
+        return f, (2 * np.cos(2 * x)) ** 2 * f - 4 * np.sin(2 * x) * f
+
+    errs = []
+    for n in (41, 81, 161):
+        x = _graded_grid(n)
+        f, exact = f_and_deriv(x)
+        op = Diff(0, x, acc=acc) ** deriv if deriv > 1 else Diff(0, x, acc=acc)
+        got = op(f)
+        b = deriv + acc + 1  # skip a few boundary points
+        errs.append(np.max(np.abs(got[b:-b] - exact[b:-b])))
+    orders = np.log2(np.array(errs[:-1]) / np.array(errs[1:]))
+    assert orders.min() > acc - 0.4
+
+
+def test_nonuniform_even_derivative_interior_matches_boundary_order():
+    """Interior even-derivative error is no worse than the boundary error.
+
+    The undersized central stencil made the interior *less* accurate than the
+    one-sided boundary stencils, which is backwards.
+    """
+    x = _graded_grid(60)
+    f = np.exp(np.sin(2 * x))
+    exact = (2 * np.cos(2 * x)) ** 2 * f - 4 * np.sin(2 * x) * f
+    got = (Diff(0, x, acc=2) ** 2)(f)
+    err = np.abs(got - exact)
+    interior = err[5:-5].max()
+    boundary = max(err[:5].max(), err[-5:].max())
+    assert interior <= boundary

tests/test_nonuniform_optimized.py

@@ -101,15 +101,17 @@ def test_batched_result_structure(self):
         assert result["center"]["coefficients"].shape[0] == len(x) - 2 * num_bndry
 
     def test_batched_even_derivative(self):
-        """Even derivatives have forward/backward stencils wider than central."""
+        """Even-derivative central stencils get no symmetry bonus on a
+        non-uniform grid, so they carry deriv+acc points (one more than the
+        one-sided schemes) to reach order acc."""
         x = np.linspace(0, 5, 20)
         result = calc_coefs_non_uni_batched(deriv=2, acc=2, coords=x)
 
         central_size = result["center"]["coefficients"].shape[1]
         forward_size = result["forward"]["coefficients"].shape[1]
 
-        # For even derivatives, forward/backward have one more point
-        assert forward_size == central_size + 1
+        assert central_size == forward_size + 1
+        assert central_size >= 2 + 2  # deriv + acc points needed for order acc
 
     def test_batched_invalid_inputs(self):
         """Invalid deriv/acc raise ValueError."""

tests/test_vector.py

@@ -196,13 +196,24 @@ def test_1d_non_uniform(self):
         laplace = Laplacian(coords=[x], acc=4)
         assert_array_almost_equal(laplace(f), expected, decimal=4)
 
-    def test_matches_uniform_result(self):
-        """Non-uniform Laplacian on a uniform grid should match the uniform version."""
+    def test_uniform_grid_gives_correct_laplacian(self):
+        """On a uniform grid the coords= Laplacian is correct and at least as
+        accurate as the h= path. Its central stencils are wider (the non-uniform
+        scheme cannot rely on the symmetry bonus), so the two paths agree only to
+        the uniform scheme's order, not bit-for-bit."""
         axes, h, [X, Y] = init_mesh(2, (50, 50))
         f = np.sin(X) * np.cos(Y)
+        expected = -2 * np.sin(X) * np.cos(Y)
         result_uniform = Laplacian(h=h, acc=2)(f)
         result_nonuniform = Laplacian(coords=axes, acc=2)(f)
-        assert_array_almost_equal(result_nonuniform, result_uniform)
+
+        interior = (slice(2, -2), slice(2, -2))
+        assert_array_almost_equal(
+            result_nonuniform[interior], expected[interior], decimal=3
+        )
+        err_uniform = np.max(np.abs(result_uniform[interior] - expected[interior]))
+        err_nonuniform = np.max(np.abs(result_nonuniform[interior] - expected[interior]))
+        assert err_nonuniform <= err_uniform
 
     def test_h_and_coords_raises(self):
         axes, h, _ = init_mesh(2, (10, 10))