TL: added derivation matrix for LagrangeApproximation + bump version to 0.1.14

Author: tlunet

CITATION.cff

@@ -16,9 +16,9 @@ authors:
     orcid: https://orcid.org/0000-0002-3879-1210
     affiliation: "Jülich Supercomputing Centre, Forschungszentrum Jülich GmbH, 52425 Jülich, Germany"
 
-version: 0.1.13
+version: 0.1.14
 doi: 10.5281/zenodo.11956479
-date-released: 2024-11-13
+date-released: 2024-11-29
 keywords: 
   - "time integration"
   - "spectral deferred corrections"

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "qmat"
-version = "0.1.13"
+version = "0.1.14"
 description = "Generation of Q-coefficients for Spectral Deferred Corrections (and other time-integration methods ...)"
 dependencies = [
     "numpy",

qmat/lagrange.py

@@ -324,3 +324,40 @@ def getIntegrationMatrix(self, intervals, numQuad='FEJER'):
         PInter = integrand.sum(axis=-1).T
 
         return PInter
+
+    def getDerivationMatrix(self):
+        r"""
+        Generate the first order derivation matrix :math:`D^{(1)}` based on the
+        Lagrange interpolant, such that
+
+        .. math::
+            D^{(1)} u \simeq \frac{du}{dx}
+
+        on the interpolation points. The formula is :
+
+        .. math::
+            D^{(1)}_{ij} = \frac{w_j/w_i}{x_i-x_j}
+
+        for :math:`i \neq j` and
+
+        .. math::
+            D^{(1)}_{jj} = -\sum_{i \neq j} D^{(1)}_{ij}`
+
+        Returns
+        -------
+        D1 : np.2darray
+            Derivation matrix.
+        """
+        w = self.weights
+        x = self.points
+
+        with np.errstate(divide='ignore'):
+            iDiff = 1 / (x[:, None] - x[None, :])
+        iDiff[np.isinf(iDiff)] = 0
+
+        base = w[None, :]/w[:, None]
+        base *= iDiff
+
+        D1 = base
+        np.fill_diagonal(D1, -D1.sum(axis=-1))
+        return D1

tests/test_2_lagrange.py

@@ -49,7 +49,7 @@ def testAsymptoticWeightsSTABLE(pType, scaleRef):
 def testInterpolation(nNodes, weightComputation):
     nodes = np.sort(np.random.rand(nNodes))
     approx = LagrangeApproximation(nodes, weightComputation=weightComputation)
-    
+
     times = np.random.rand(nNodes*2)
     P = approx.getInterpolationMatrix(times)
 
@@ -66,7 +66,7 @@ def testEvaluation(nNodes, weightComputation):
     times = np.random.rand(nNodes*2)
     polyCoeffs = np.random.rand(nNodes)
     polyValues = np.polyval(polyCoeffs, nodes)
-    
+
     approx = LagrangeApproximation(nodes, weightComputation=weightComputation, fValues=polyValues)
     P = approx.getInterpolationMatrix(times)
     refEvals = P @ polyValues
@@ -76,7 +76,7 @@ def testEvaluation(nNodes, weightComputation):
 
     polyEvals = approx(t=times)
     assert np.allclose(polyEvals, refEvals)
-    
+
 
 @pytest.mark.parametrize("numQuad", ["LEGENDRE_NUMPY", "LEGENDRE_SCIPY", "FEJER"])
 @pytest.mark.parametrize("weightComputation", ["AUTO", "FAST", "STABLE", "CHEBFUN"])
@@ -92,4 +92,19 @@ def testIntegration(nNodes, weightComputation, numQuad):
     polyNodes = np.polyval(polyCoeffs, nodes)
     polyInteg = np.polyval(np.polyint(polyCoeffs), times) - np.polyval(np.polyint(polyCoeffs), 0)
 
-    assert np.allclose(polyInteg, P @ polyNodes)
+    assert np.allclose(polyInteg, P @ polyNodes)
+
+
+@pytest.mark.parametrize("weightComputation", ["AUTO", "FAST", "STABLE", "CHEBFUN"])
+@pytest.mark.parametrize("nNodes", nNodeTests)
+def testDerivation(nNodes, weightComputation):
+    nodes = np.sort(np.random.rand(nNodes))
+    approx = LagrangeApproximation(nodes, weightComputation=weightComputation)
+
+    D = approx.getDerivationMatrix()
+
+    polyCoeffs = np.random.rand(nNodes)
+    polyNodes = np.polyval(polyCoeffs, nodes)
+    polyDeriv = np.polyval(np.polyder(polyCoeffs), nodes)
+
+    assert np.allclose(polyDeriv, D @ polyNodes)