Look at - https://galgebra.readthedocs.io/en/latest/
gradient (geometric derivative) and associated operators implemented with the same interface. Nothing to do with integration implemented. Attached is the code for Dop.py and the test code in Jupyter notebook. The problem with the differential operator class is that in order to use them to stuff a sympy Matrix they have to be sympy expressions and I don't know if that is possible without changing some of the sympy core. You might want to look at geometric calculus -
https://en.wikipedia.org/wiki/Geometric_calculusgalgebra includes all of the differential geometric calculus. Could you give me examples of what you want to do in pdf/latex format files. Essentially all that Dop.py does is provide a different api to using sympy derivatives. If all is required is changing the api that would not be too difficult. A integrator class could be written with a more natural api. What would be difficult is if you had an integrand with unevaluated differential expressions. I don't think sympy can return f for the integral of (df/dx)dx without first differentiating and then integrating.
On 4/19/22 12:54 PM, Jonathan Gutow wrote:
Alan,I have thought about this a little too. I have not had time to work on it recently. The issue I ran into is that to make this work well in SymPy you really need the concept of an infinitesimal dx, dy, dz, etc. Things got circular when I tried to implement that using the sympy definition of functions. It might be worth comparing notes. I would really like to be able to work with total differentials and do integrations on expressions containing infinitesimals in the Leibniz style notation.JonathanOn Apr 19, 2022, at 10:24 AM, Alan Bromborsky <[email protected]> wrote:CAUTION:This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.I am working on a linear differential operator class for sympy (Dop.py). Here is a Jupyter notebook with a simple example -<sdop_test.jpg> Note that lap*f gives a function but f*lap gives differential operator.--You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/7C0243F5-5603-4DBF-806D-BC2C51507E86%40uwosh.edu <https://groups.google.com/d/msgid/sympy/7C0243F5-5603-4DBF-806D-BC2C51507E86%40uwosh.edu?utm_medium=email&utm_source=footer>.
-- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/31039ccc-0ca6-ca30-d5d0-9f36ec6107a4%40gmail.com.
"""
Differential operators, for all sympy expressions
For multivector-customized differential operators, see :class:`galgebra.mv.Dop`.
"""
import copy
import itertools
import numbers
import warnings
from typing import List, Tuple, Any, Iterable
from sympy import Symbol, S, Add, simplify, diff, Expr, Dummy, Matrix, shape
from IPython.display import display, Math, Latex
ZERO_STR = ' 0 '
PD = '\u2202'
PDpow = ('','','\u00B2','\u00B3','\u2074','\u2075','\u2076','\u2077','\u2078','\u2079')
PDmax = 9
def apply_function_list(f, x):
if isinstance(f, (tuple, list)):
fx = x
for fi in f:
fx = fi(fx)
return fx
else:
return f(x)
class Simp:
modes = [simplify]
@staticmethod
def profile(s):
Simp.modes = s
@staticmethod
def apply(expr):
obj = S.Zero
for coef, base in linear_expand_terms(expr):
obj += apply_function_list(Simp.modes, coef) * base
return obj
@staticmethod
def applymv(mv):
return Mv(Simp.apply(mv.obj), ga=mv.Ga)
def _consolidate_terms(terms):
"""
Remove zero coefs and consolidate coefs with repeated pdiffs.
"""
new_coefs = []
new_pdiffs = []
for coef, pd in terms:
if coef != S.Zero:
if pd in new_pdiffs:
index = new_pdiffs.index(pd)
new_coefs[index] += coef
else:
new_coefs.append(coef)
new_pdiffs.append(pd)
return tuple(zip(new_coefs, new_pdiffs))
def _merge_terms(terms1, terms2):
""" Concatenate and consolidate two sets of already-consolidated terms """
pdiffs1 = [pdiff for _, pdiff in terms1]
pdiffs2 = [pdiff for _, pdiff in terms2]
pdiffs = pdiffs1 + [x for x in pdiffs2 if x not in pdiffs1]
coefs = len(pdiffs) * [S.Zero]
for coef, pdiff in terms1:
index = pdiffs.index(pdiff)
coefs[index] += coef
for coef, pdiff in terms2:
index = pdiffs.index(pdiff)
coefs[index] += coef
# remove zeros
return [(coef, pdiff) for coef, pdiff in zip(coefs, pdiffs) if coef != S.Zero]
def _eval_derivative_n_times_terms(terms, x, n):
for i in range(n):
new_terms = []
for k, term in enumerate(terms):
dc = _basic_diff(term[0], x)
pd = _basic_diff(term[1], x)
# print 'D0, term, dc, pd =', D0, term, dc, pd
if dc != 0:
new_terms.append((dc, term[1]))
if pd != 0:
new_terms.append((term[0], pd))
terms = new_terms
return _consolidate_terms(terms)
################ Scalar Partial Differential Operator Class ############
class _BaseDop():
""" Base class for differential operators - used to avoid accidental promotion """
pass
class Sdop(_BaseDop):
"""
Scalar differential operator is of the form (Einstein summation)
.. math:: D = c_{i}*D_{i}
where the :math:`c_{i}`'s are scalar coefficient (they could be functions)
and the :math:`D_{i}`'s are partial differential operators (:class:`Pdop`).
Attributes
----------
terms : tuple of tuple
the structure :math:`((c_{1},D_{1}),(c_{2},D_{2}), ...)`
"""
def TSimplify(self):
return Sdop([
(Simp.apply(coef), pdiff) for coef, pdiff in self.terms
])
@staticmethod
def consolidate_coefs(sdop):
"""
Remove zero coefs and consolidate coefs with repeated pdiffs.
"""
if isinstance(sdop, Sdop):
return Sdop(_consolidate_terms(sdop.terms))
else:
return _consolidate_terms(sdop)
def simplify(self, modes=simplify):
return Sdop([
(apply_function_list(modes, coef), pdiff)
for coef, pdiff in self.terms
])
def _with_sorted_terms(self):
new_terms = sorted(self.terms, key=lambda term: Pdop.sort_key(term[1]))
return Sdop(new_terms)
def __str__(self):
if len(self.terms) == 0:
return ZERO_STR
self = self._with_sorted_terms()
s = ''
for coef, pdop in self.terms:
coef_str = str(coef)
pd_str = str(pdop)
if coef == S.One:
s += pd_str
elif coef == S.NegativeOne:
s += '-' + pd_str
else:
if isinstance(coef, Add):
s += '(' + coef_str + ')*' + pd_str
else:
s += coef_str + '*' + pd_str
s += ' + '
s = s.replace('+ -', '- ')
s = s[:-3]
return s
def _repr_latex_(self):
if len(self.terms) == 0:
return ZERO_STR
self = self._with_sorted_terms()
s = ''
for coef, pdop in self.terms:
coef_str = str(coef)
pd_str = pdop._repr_latex_()
if coef == S.One:
if pd_str == '':
s += '1'
else:
s += pd_str
elif coef == S.NegativeOne:
if pd_str == '':
s += '-1'
else:
s += '-' + pd_str
else:
if isinstance(coef, Add):
s += r'\left ( ' + coef_str + r'\right ) ' + pd_str
else:
s += coef_str + ' ' + pd_str
s += ' + '
s = s.replace('+ -', '- ')
return s[:-3]
def __init_from_symbol(self, symbol: Symbol) -> None:
self.terms = ((S.One, Pdop(symbol)),)
def __init_from_coef_and_pdiffs(self, coefs: List[Any], pdiffs: List['Pdop']) -> None:
if not isinstance(coefs, list) or not isinstance(pdiffs, list):
raise TypeError("coefs and pdiffs must be lists")
if len(coefs) != len(pdiffs):
raise ValueError('In Sdop.__init__ coefficent list and Pdop list must be same length.')
self.terms = tuple(zip(coefs, pdiffs))
def __init_from_terms(self, terms: Iterable[Tuple[Any, 'Pdop']]) -> None:
self.terms = tuple(terms)
def __init__(self, *args):
if len(args) == 1:
if isinstance(args[0], Symbol):
self.__init_from_symbol(*args)
elif isinstance(args[0], (list, tuple)):
self.__init_from_terms(*args)
else:
raise TypeError(
"A symbol or sequence is required (got type {})"
.format(type(args[0]).__name__))
elif len(args) == 2:
self.__init_from_coef_and_pdiffs(*args)
else:
raise TypeError(
"Sdop() takes from 1 to 2 positional arguments but {} were "
"given".format(len(args)))
def __call__(self, arg):
# Ensure that we return the right type even when there are no terms - we
# do this by adding `0 * d(arg)/d(nonexistant)`, which must be zero, but
# will be a zero of the right type.
dummy_var = Dummy('nonexistant')
terms = self.terms or ((S.Zero, Pdop(dummy_var)),)
return sum([coef * pdiff(arg) for coef, pdiff in terms])
def __neg__(self):
return Sdop([(-coef, pdiff) for coef, pdiff in self.terms])
@staticmethod
def Add(sdop1, sdop2):
if isinstance(sdop1, Sdop) and isinstance(sdop2, Sdop):
return Sdop(_merge_terms(sdop1.terms, sdop2.terms))
else:
# convert values to multiplicative operators
if not isinstance(sdop2, _BaseDop):
sdop2 = Sdop([(sdop2, Pdop({}))])
elif not isinstance(sdop1, _BaseDop):
sdop1 = Sdop([(sdop1, Pdop({}))])
else:
return NotImplemented
return Sdop.Add(sdop1, sdop2)
def __eq__(self, other):
if isinstance(other, Sdop):
diff = self - other
return len(diff.terms) == 0
else:
return NotImplemented
def __add__(self, sdop):
return Sdop.Add(self, sdop)
def __radd__(self, sdop):
return Sdop.Add(sdop, self)
def __sub__(self, sdop):
return Sdop.Add(self, -sdop)
def __rsub__(self, sdop):
return Sdop.Add(-self, sdop)
def __mul__(self, sdopr):
# alias for applying the operator
return self.__call__(sdopr)
def __rmul__(self, sdop):
return Sdop([(sdop * coef, pdiff) for coef, pdiff in self.terms])
def _eval_derivative_n_times(self, x, n):
return Sdop(_eval_derivative_n_times_terms(self.terms, x, n))
def __pow__(self, n):
if ( isinstance(n,int)):
if ( n == 1 ):
return self
x = self
for _ in itertools.repeat(None, n-1):
x = self*x
return x
else:
return NotImplemented
#################### Partial Derivative Operator Class #################
def _basic_diff(f, x, n=1):
""" Simple wrapper for `diff` that works for our types too """
if isinstance(f, (Expr, Symbol, numbers.Number)): # f is sympy expression
return diff(f, x, n)
elif hasattr(f, '_eval_derivative_n_times'):
# one of our types
return f._eval_derivative_n_times(x, n)
else:
raise ValueError('In_basic_diff type(arg) = ' + str(type(f)) + ' not allowed.')
class Pdop(_BaseDop):
r"""
Partial derivative operatorp.
The partial derivatives are of the form
.. math::
\partial_{i_{1}...i_{n}} =
\frac{\partial^{i_{1}+...+i_{n}}}{\partial{x_{1}^{i_{1}}}...\partial{x_{n}^{i_{n}}}}.
If :math:`i_{j} = 0` then the partial derivative does not contain the
:math:`x^{i_{j}}` coordinate.
Attributes
----------
pdiffs : dict
A dictionary where coordinates are keys and key value are the number of
times one differentiates with respect to the key.
order : int
Total number of differentiations.
When this is zero (i.e. when :attr:`pdiffs` is ``{}``) then this object
is the identity operator, and returns its operand unchanged.
"""
def sort_key(self, order=None):
return (
# lower order derivatives first
self.order,
# sorted by symbol after that, after expansion
sorted([
x.sort_key(order)
for x, k in self.pdiffs.items()
for i in range(k)
])
)
def __eq__(self, A):
if isinstance(A, Pdop) and self.pdiffs == A.pdiffs:
return True
else:
if len(self.pdiffs) == 0 and A == S.One:
return True
return False
def __init__(self, __arg):
"""
The partial differential operator is a partial derivative with
respect to a set of real symbols (variables).
"""
# galgebra 0.4.5
if __arg is None:
warnings.warn(
"`Pdop(None)` is deprecated, use `Pdop({})` instead",
DeprecationWarning, stacklevel=2)
__arg = {}
if isinstance(__arg, dict): # Pdop defined by dictionary
self.pdiffs = __arg
elif isinstance(__arg, Symbol): # First order derivative with respect to symbol
self.pdiffs = {__arg: 1}
else:
raise TypeError('A dictionary or symbol is required, got {!r}'.format(__arg))
self.order = sum(self.pdiffs.values())
def _eval_derivative_n_times(self, x, n) -> 'Pdop': # pdiff(self)
# d is partial derivative
pdiffs = copy.copy(self.pdiffs)
if x in pdiffs:
pdiffs[x] += n
else:
pdiffs[x] = n
return Pdop(pdiffs)
def __call__(self, arg):
"""
Calculate nth order partial derivative (order defined by
self) of expression
"""
for x, n in self.pdiffs.items():
arg = _basic_diff(arg, x, n)
return arg
def __mul__(self, other): # functional product of self and arg (self*arg)
return self(other)
def __rmul__(self, other): # functional product of arg and self (arg*self)
assert not isinstance(other, Pdop)
return Sdop([(other, self)])
def __str__(self):
if self.order == 0:
return '1'
s = ''
for x in self.pdiffs:
s += '('
n = self.pdiffs[x]
if n == 1:
s += PD
else:
s += PD+'^'+str(n)
s += '/' + PD + str(x) + ')'
return s
def _repr_latex_(self):
if self.order == 0:
return ''
s = r'\frac{\partial'
if self.order > 1:
s += '^{' + str(self.order) + '}'
s += '}{'
keys = list(self.pdiffs.keys())
keys.sort(key=lambda x: x.sort_key())
for key in keys:
i = self.pdiffs[key]
s += r'\partial ' + str(key)
if i > 1:
s += '^{' + str(i) + '}'
s += '}'
return s
def gprint(*args,strng = False):
s = ''
lp = r' \left ( '
rp = r' \right ) '
lb = r' \left [ '
rb = r' \right ] '
lft = ''
rght = ''
for arg in args:
if isinstance(arg, str):
s += arg
elif isinstance(arg, (tuple,list)):
if isinstance(arg, tuple):
lft = lp
rght = rp
else:
lft = lb
rght = rb
s += lft
for iarg in arg:
s+= gprint(iarg,strng=True)+','
s = s[:-1] + rght
elif isinstance(arg, Matrix):
ncol = shape(arg)[1]
nrow = shape(arg)[0]
s += lb+r'\begin{array}{'+ncol*'c'+'}'
for irow in range(nrow):
for icol in range(ncol):
s += arg[irow,icol]._repr_latex_() + ' & '
s = s[:-3]+r'\\'
s = s[:-2]+r'\end{array}'+rb
else:
s += arg._repr_latex_()
if strng:
return s
else:
s = s.replace('$','')
display(Math(s))
return
sdop_test.ipynb
Description: application/ipynb
