Merge pull request #74 from maroba/symbolic

Version 0.10.0: Introduce symbolic capabilities.

README.md

@@ -18,6 +18,7 @@ any number of dimensions.
 * Fully vectorized for speed
 * Matrix representations of arbitrary linear differential operators
 * Solve partial differential equations with Dirichlet or Neumann boundary conditions
+* **New in version 0.10**: Symbolic representation of finite difference schemes
 
 ## Installation
 
@@ -216,6 +217,30 @@ which returns
 {(0, 0): -2.0, (1, 1): 0.5, (-1, -1): 0.5, (1, -1): 0.5, (-1, 1): 0.5}
 ```
 
+## Symbolic Representations
+
+As of version 0.10, findiff can also provide a symbolic representation of finite difference schemes suitable for using in conjunction with sympy. The main use case is to facilitate deriving your own iteration schemes.
+
+```python
+from findiff import SymbolicMesh, SymbolicDiff
+
+mesh = SymbolicMesh("x, y")
+d2_dx2, d2_dy2 = [SymbolicDiff(mesh, axis=k, degree=2) for k in range(2)]
+
+(
+    d2_dx2(u, at=(m, n), offsets=(-1, 0, 1)) + 
+    d2_dy2(u, at=(m, n), offsets=(-1, 0, 1))
+)
+```
+
+Outputs:
+
+$$
+\frac{{u}_{m,n + 1} + {u}_{m,n - 1} - 2 {u}_{m,n}}{\Delta y^{2}} + \frac{{u}_{m + 1,n} + {u}_{m - 1,n} - 2 {u}_{m,n}}{\Delta x^{2}}
+$$
+
+Also see the [example notebook](examples/symbolic.ipynb).
+
 ## Partial Differential Equations
 
 _findiff_ can be used to easily formulate and solve partial differential equation problems

findiff/__init__.py

@@ -14,11 +14,13 @@
 - New in version 0.7: Generate matrix representations of arbitrary linear differential operators
 - New in version 0.8: Solve partial differential equations with Dirichlet or Neumann boundary conditions
 - New in version 0.9: Generate differential operators for generic stencils
+- New in version 0.10: Create symbolic representations of finite difference schemes
 """
 
-__version__ = '0.9.2'
+__version__ = '0.10.0'
 
 from .coefs import coefficients
 from .operators import FinDiff, Coef, Identity, Coefficient
 from .vector import Gradient, Divergence, Curl, Laplacian
 from .pde import PDE, BoundaryConditions
+from .symbolic import SymbolicMesh, SymbolicDiff

findiff/symbolic.py

@@ -0,0 +1,140 @@
+from sympy import IndexedBase, Symbol, Add, Indexed, latex
+from sympy import Rational, factorial, linsolve, Matrix
+
+
+class SymbolicMesh:
+    """
+    Represents the mesh on which to evaluate finite difference approximations.
+    """
+
+    def __init__(self, coord, equidistant=True):
+        """Constructor.
+
+        Parameters
+        ----------
+        coord: str
+            A comma-separated string of coordinate names for the mesh,
+            e.g. "x,y" or simply "x"
+
+        equidistant: bool
+            Flag indicating whether the mesh is equidistant.
+        """
+        assert isinstance(coord, str)
+
+        self.equidistant = equidistant
+        self._coord_names = [n.replace(" ", "") for n in coord.split(",")]
+        self._coord = [IndexedBase(name) for name in self._coord_names]
+
+    @property
+    def ndims(self):
+        """The number of dimensions of the mesh."""
+        return len(self._coord)
+
+    @property
+    def coord(self):
+        """
+        Returns a tuple with the symbols for the coordinates.
+        """
+        return self._coord
+
+    @property
+    def spacing(self):
+        """
+        Returns a tuple with the spacing of the mesh along all axes.
+        Only makes sense for equidistant grid. Raises exception in
+        case of non-equidistant grids.
+        """
+        if self.equidistant:
+            spacings = tuple(Symbol(f"\\Delta {x}") for x in self._coord_names)
+            return spacings
+        raise Exception("Non-equidistant mesh does not have spacing property.")
+
+    @staticmethod
+    def create_symbol(name):
+        """
+        Creates a *sympy* symbol of a given name which can carry as many
+        indices as the mesh has dimensions.
+
+        Parameters
+        ----------
+        name: str
+            The name of the meshed symbol.
+
+        Returns
+        -------
+        An index-carrying *sympy* symbol (IndexedBase).
+        """
+        return IndexedBase(name)
+
+
+class SymbolicDiff:
+    """
+    A symbolic representation of the finite difference approximation
+    of a partial derivative. Based on *sympy*.
+    """
+
+    def __init__(self, mesh, axis=0, degree=1):
+        """Constructor
+
+        Parameters
+        ----------
+        mesh: SymbolicMesh
+            The symbolic grid on which to evaluate the derivative.
+        axis: int
+            The index of the axis with respect to which to differentiate.
+        degree: int > 0
+            The degree of the partial derivative.
+
+        """
+        self.mesh = mesh
+        self.axis = axis
+        self.degree = degree
+
+    def __call__(self, f, at, offsets):
+        if not isinstance(at, tuple) and not isinstance(at, list):
+            at = [at]
+
+        if self.mesh.ndims != len(at):
+            raise ValueError("Index tuple must match the number of dimensions!")
+
+        coefs = self._compute_coefficients(f, at, offsets)
+        terms = []
+        for coef, off in zip(coefs, offsets):
+            inds = list(at)
+            inds[self.axis] += off
+            inds = tuple(inds)
+            terms.append(coef * f[inds])
+
+        return Add(*terms).simplify()
+
+    def _compute_coefficients(self, f, at, offsets):
+
+        n = len(offsets)
+        # the first row always contains 1s:
+        matrix = [[1] * n]
+
+        def spac(off):
+            """A helper function to get the spacing between grid points."""
+            if self.mesh.equidistant:
+                h = self.mesh.spacing[self.axis]
+            else:
+                x = self.mesh.coord[self.axis]
+                h = x[at[self.axis] + off] - x[at[self.axis]]
+            return h
+
+        # build up the matrix incrementally:
+        for i in range(1, n):
+            ifac = Rational(1, factorial(i))
+            row = [ifac * (off * spac(off)) ** i for off in offsets]
+            matrix.append(row)
+
+        # only the entry corresponding to the requested derivative degree
+        # is non-zero:
+        rhs = [0] * n
+        rhs[self.degree] = 1
+
+        # solve the equation system
+        matrix = Matrix(matrix)
+        rhs = Matrix(rhs)
+        sol = linsolve((matrix, rhs))
+        return list(sol)[0].simplify()

setup.py

@@ -39,7 +39,7 @@ def get_version():
         * _New in version 0.7:_ Generate matrix representations of arbitrary linear differential operators
         * _New in version 0.8:_ Solve partial differential equations with Dirichlet or Neumann boundary conditions
         * _New in version 0.9:_ Generate differential operators for generic stencils
-
+        * _New in version 0.10:_ Create symbolic representations of finite difference schemes
     """,
 
     license='MIT',

test/test_symbolic.py

@@ -0,0 +1,101 @@
+import unittest
+
+from sympy import IndexedBase, Symbol, Expr, Eq, symbols, latex
+
+from findiff.symbolic import SymbolicMesh, SymbolicDiff
+
+
+class TestSymbolicMesh(unittest.TestCase):
+
+    def test_parse_symbolic_mesh(self):
+        # 1D
+        mesh = SymbolicMesh(coord="x", equidistant=True)
+        x, = mesh.coord
+        dx, = mesh.spacing
+
+        self.assertEqual(IndexedBase, type(x))
+        self.assertEqual(Symbol, type(dx))
+
+        # 2D
+        mesh = SymbolicMesh(coord="x,y", equidistant=True)
+        x, y = mesh.coord
+        dx, dy = mesh.spacing
+
+        self.assertEqual(IndexedBase, type(x))
+        self.assertEqual(IndexedBase, type(y))
+        self.assertEqual(Symbol, type(dx))
+        self.assertEqual(Symbol, type(dy))
+
+        # ignores whitespace
+        mesh = SymbolicMesh(coord="x,  y", equidistant=True)
+        x, y = mesh.coord
+
+        self.assertEqual("x", str(x))
+        self.assertEqual("y", str(y))
+
+    def test_create_symbol(self):
+        # defaults
+        mesh = SymbolicMesh(coord="x", equidistant=True)
+        actual = mesh.create_symbol("u")
+        expected = IndexedBase("u")
+        self.assertEqual(latex(actual), latex(expected))
+
+        # both indices down
+        mesh = SymbolicMesh(coord="x,y", equidistant=True)
+        n, m = symbols("n, m")
+        actual = latex(mesh.create_symbol("u")[n, m])
+        #expected = "u{}_{n}{}_{m}"
+        expected = "{u}_{n,m}"
+        self.assertEqual(latex(actual), latex(expected))
+
+        # both indices up
+        #mesh = SymbolicMesh(coord="x,y", equidistant=True)
+        #n, m = symbols("n, m")
+        #u = mesh.create_symbol("u", pos="u,u")
+        #actual = latex(u[n, m])
+        #expected = "u{}^{n}{}^{m}"
+        #self.assertEqual(actual, expected)
+
+
+class TestDiff(unittest.TestCase):
+
+    def test_init(self):
+        mesh = SymbolicMesh("x")
+        d = SymbolicDiff(mesh)
+
+        self.assertEqual(d.axis, 0)
+        self.assertEqual(d.degree, 1)
+        self.assertEqual(id(mesh), id(d.mesh))
+
+    def test_call(self):
+        # 1D
+        mesh = SymbolicMesh("x")
+        u = mesh.create_symbol("u")
+        d = SymbolicDiff(mesh)
+        n = Symbol("n")
+
+        actual = d(u, at=n, offsets=[-1, 0, 1])
+
+        expected = (u[n + 1] - u[n - 1]) / (2 * mesh.spacing[0])
+
+        self.assertEqual(
+            0, (expected - actual).simplify()
+        )
+
+        # 2D
+        mesh = SymbolicMesh("x, y")
+        u = mesh.create_symbol("u")
+        d = SymbolicDiff(mesh, axis=1)
+        n, m = symbols("n, m")
+
+        actual = d(u, at=(n, m), offsets=[-1, 0, 1])
+
+        expected = (u[n, m + 1] - u[n, m - 1]) / (2 * mesh.spacing[1])
+
+        self.assertEqual(
+            0, (expected - actual).simplify()
+        )
+
+
+if __name__ == '__main__':
+    unittest.main()