Rm findiff

pyproject.toml

@@ -35,7 +35,6 @@ classifiers = [
 ]
 requires-python = ">=3.10"
 dependencies = [
-    "findiff>=0.9.2",
     "h5py>=3.8.0",
     "numpy>=1.21.5",
     "scipy>=1.9.0",

src/WallGo/boltzmann.py

@@ -8,14 +8,14 @@
 import logging
 from enum import Enum, auto
 import numpy as np
-import findiff  # finite difference methods
 from .containers import BoltzmannBackground, BoltzmannDeltas
 from .grid import Grid
 from .polynomial import Polynomial, SpectralConvergenceInfo
 from .particle import Particle
 from .collisionArray import CollisionArray
 from .results import BoltzmannResults
 from .exceptions import CollisionLoadError
+from .helpers import finiteDiffMatrix
 
 if typing.TYPE_CHECKING:
     import importlib
@@ -704,10 +704,8 @@ def buildLinearEquations(
             intertwinerRpMat = np.identity(self.grid.N - 1)
             # derivative matrices
             chiFull, rzFull, _ = self.grid.getCompactCoordinates(endpoints=True)
-            derivOperatorChi = findiff.FinDiff((0, chiFull, 1), acc=2)
-            derivMatrixChi = derivOperatorChi.matrix((self.grid.M + 1,))
-            derivOperatorRz = findiff.FinDiff((0, rzFull, 1), acc=2)
-            derivMatrixRz = derivOperatorRz.matrix((self.grid.N + 1,))
+            derivMatrixChi = finiteDiffMatrix(chiFull)
+            derivMatrixRz = finiteDiffMatrix(rzFull)
             # spatial derivatives of profiles, endpoints used for taking
             # derivatives but then dropped as deltaF fixed at 0 at endpoints
             dTemperaturedChi = (derivMatrixChi @ temperatureFull)[

src/WallGo/helpers.py

@@ -425,6 +425,34 @@ def hessian(
     fEvaluation = f(pos, *args).reshape(shape)
     return np.asarray(np.sum(coeff * fEvaluation, axis=-3))
 
+def finiteDiffMatrix(grid: np.ndarray) -> np.ndarray:
+    """
+    Computes the finite difference matrix corresponding to the array grid.
+    grid must contain the endpoints.
+    """
+    diffMatrix = np.zeros((grid.size, grid.size))
+    dx1 = -(grid[1:-1] - grid[:-2])
+    dx2 = grid[2:] - grid[1:-1]
+
+    # Below the diagonal
+    diffMatrix[1:-1,:-2] += np.diag(dx2/(dx1*(dx2-dx1)))
+    # Diagonal
+    diffMatrix[1:-1,1:-1] += np.diag(-(dx1+dx2)/(dx1*dx2))
+    # Above the diagonal
+    diffMatrix[1:-1,2:] += np.diag(dx1/(dx2*(dx1-dx2)))
+
+    # Left boundary
+    diffMatrix[0,0] = ((2*grid[0] - grid[1] - grid[2])/((grid[0] - grid[1])*(grid[0] - grid[2])))
+    diffMatrix[0,1] = ((-grid[0] + grid[2])/((grid[0] - grid[1])*(grid[1] - grid[2])))
+    diffMatrix[0,2] = ((-grid[0] + grid[1])/((grid[0] - grid[2])*(-grid[1] + grid[2])))
+
+    # Right boundary
+    diffMatrix[-1,-1] = ((2*grid[-1] - grid[-2] - grid[-3])/((grid[-1] - grid[-2])*(grid[-1] - grid[-3])))
+    diffMatrix[-1,-2] = ((-grid[-1] + grid[-3])/((grid[-1] - grid[-2])*(grid[-2] - grid[-3])))
+    diffMatrix[-1,-3] = ((-grid[-1] + grid[-2])/((grid[-1] - grid[-3])*(-grid[-2] + grid[-3])))
+
+    return diffMatrix
+
 
 def gammaSq(v: float) -> float:
     r"""