From 6259d04d11ce6fd516a4c504399064da5e72e98a Mon Sep 17 00:00:00 2001 From: maroba Date: Tue, 18 Jul 2023 16:50:48 +0200 Subject: [PATCH] Version 0.10.0 --- README.md | 25 ++++++++ findiff/__init__.py | 4 +- findiff/symbolic.py | 140 ++++++++++++++++++++++++++++++++++++++++++ setup.py | 2 +- test/test_symbolic.py | 101 ++++++++++++++++++++++++++++++ 5 files changed, 270 insertions(+), 2 deletions(-) create mode 100644 findiff/symbolic.py create mode 100644 test/test_symbolic.py diff --git a/README.md b/README.md index 41ac88d..d67878d 100644 --- a/README.md +++ b/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 diff --git a/findiff/__init__.py b/findiff/__init__.py index 3fe5a64..d6ada2c 100644 --- a/findiff/__init__.py +++ b/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 diff --git a/findiff/symbolic.py b/findiff/symbolic.py new file mode 100644 index 0000000..b90f676 --- /dev/null +++ b/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() diff --git a/setup.py b/setup.py index 36b933f..ee7cb71 100644 --- a/setup.py +++ b/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', diff --git a/test/test_symbolic.py b/test/test_symbolic.py new file mode 100644 index 0000000..83300c2 --- /dev/null +++ b/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()