Operador diferencial utilizable en forma de matriz, en el módulo Sympy de Python

Necesitamos dos matrices de operadores diferenciales [B] y [C] tales como:

 B = sympy.Matrix([[ D(x), D(y) ], [ D(y), D(x) ]]) C = sympy.Matrix([[ D(x), D(y) ]]) ans = B * sympy.Matrix([[x*y**2], [x**2*y]]) print ans [x**2 + y**2] [ 4*x*y] ans2 = ans * C print ans2 [2*x, 2*y] [4*y, 4*x] 

Esto también podría aplicarse para calcular el rizo de un campo vectorial como:

 culr = sympy.Matrix([[ D(x), D(y), D(z) ]]) field = sympy.Matrix([[ x**2*y, x*y*z, -x**2*y**2 ]]) 

Para resolver esto utilizando Sympy, se tuvo que crear la siguiente clase de Python:

 import sympy class D( sympy.Derivative ): def __init__( self, var ): super( D, self ).__init__() self.var = var def __mul__(self, other): return sympy.diff( other, self.var ) 

Esta clase por sí sola se resuelve cuando la matriz de operadores diferenciales se multiplica por la izquierda. Aquí se ejecuta diff solo cuando se conoce la función a diferenciar.

Para solucionar el problema cuando la matriz de operadores diferenciales se multiplica a la derecha, el método __mul__ en la clase principal Expr tuvo que cambiarse de la siguiente manera:

 class Expr(Basic, EvalfMixin): # ... def __mul__(self, other): import sympy if other.__class__.__name__ == 'D': return sympy.diff( self, other.var ) else: return Mul(self, other) #... 

Funciona bastante bien, pero debería haber una mejor solución nativa en Sympy para manejar esto. ¿Alguien sabe lo que podría ser?

Esta solución aplica los consejos de las otras respuestas y de aquí . El operador D se puede definir de la siguiente manera:

  • solo se considera cuando se multiplica por la izquierda, de modo que D(t)*2*t**3 = 6*t**2 pero 2*t**3*D(t) no hace nada
  • todas las expresiones y símbolos utilizados con D deben tener is_commutative = False
  • se evalúa en el contexto de una expresión dada usando evaluateExpr()
    • que va de la derecha a la izquierda a lo largo de la expresión que encuentra a los operadores D y aplica mydiff() * a la parte derecha correspondiente

*: mydiff se utiliza en lugar de diff para permitir que se cree un orden D superior, como mydiff(D(t), t) = D(t,t)

La diff dentro de __mul__() en D se mantuvo solo como referencia, ya que en la solución actual, evaluateExpr() realiza el trabajo de diferenciación. Un mudule de python fue creado y guardado como d.py

 import sympy from sympy.core.decorators import call_highest_priority from sympy import Expr, Matrix, Mul, Add, diff from sympy.core.numbers import Zero class D(Expr): _op_priority = 11. is_commutative = False def __init__(self, *variables, **assumptions): super(D, self).__init__() self.evaluate = False self.variables = variables def __repr__(self): return 'D%s' % str(self.variables) def __str__(self): return self.__repr__() @call_highest_priority('__mul__') def __rmul__(self, other): return Mul(other, self) @call_highest_priority('__rmul__') def __mul__(self, other): if isinstance(other, D): variables = self.variables + other.variables return D(*variables) if isinstance(other, Matrix): other_copy = other.copy() for i, elem in enumerate(other): other_copy[i] = self * elem return other_copy if self.evaluate: return diff(other, *self.variables) else: return Mul(self, other) def __pow__(self, other): variables = self.variables for i in range(other-1): variables += self.variables return D(*variables) def mydiff(expr, *variables): if isinstance(expr, D): expr.variables += variables return D(*expr.variables) if isinstance(expr, Matrix): expr_copy = expr.copy() for i, elem in enumerate(expr): expr_copy[i] = diff(elem, *variables) return expr_copy return diff(expr, *variables) def evaluateMul(expr): end = 0 if expr.args: if isinstance(expr.args[-1], D): if len(expr.args[:-1])==1: cte = expr.args[0] return Zero() end = -1 for i in range(len(expr.args)-1+end, -1, -1): arg = expr.args[i] if isinstance(arg, Add): arg = evaluateAdd(arg) if isinstance(arg, Mul): arg = evaluateMul(arg) if isinstance(arg, D): left = Mul(*expr.args[:i]) right = Mul(*expr.args[i+1:]) right = mydiff(right, *arg.variables) ans = left * right return evaluateMul(ans) return expr def evaluateAdd(expr): newargs = [] for arg in expr.args: if isinstance(arg, Mul): arg = evaluateMul(arg) if isinstance(arg, Add): arg = evaluateAdd(arg) if isinstance(arg, D): arg = Zero() newargs.append(arg) return Add(*newargs) #courtesy: https://stackoverflow.com/a/48291478/1429450 def disableNonCommutivity(expr): replacements = {s: sympy.Dummy(s.name) for s in expr.free_symbols} return expr.xreplace(replacements) def evaluateExpr(expr): if isinstance(expr, Matrix): for i, elem in enumerate(expr): elem = elem.expand() expr[i] = evaluateExpr(elem) return disableNonCommutivity(expr) expr = expr.expand() if isinstance(expr, Mul): expr = evaluateMul(expr) elif isinstance(expr, Add): expr = evaluateAdd(expr) elif isinstance(expr, D): expr = Zero() return disableNonCommutivity(expr) 

Ejemplo 1: rizo de un campo vectorial. Tenga en cuenta que es importante definir las variables con commutative=False ya que su orden en Mul().args afectará los resultados, consulte esta otra pregunta .

 from d import D, evaluateExpr from sympy import Matrix sympy.var('x', commutative=False) sympy.var('y', commutative=False) sympy.var('z', commutative=False) curl = Matrix( [[ D(x), D(y), D(z) ]] ) field = Matrix( [[ x**2*y, x*y*z, -x**2*y**2 ]] ) evaluateExpr( curl.cross( field ) ) # [-x*y - 2*x**2*y, 2*x*y**2, -x**2 + y*z] 

Ejemplo 2: aproximación de Ritz típica utilizada en el análisis estructural.

 from d import D, evaluateExpr from sympy import sin, cos, Matrix sin.is_commutative = False cos.is_commutative = False g1 = [] g2 = [] g3 = [] sympy.var('x', commutative=False) sympy.var('t', commutative=False) sympy.var('r', commutative=False) sympy.var('A', commutative=False) m=5 n=5 for j in xrange(1,n+1): for i in xrange(1,m+1): g1 += [sin(i*x)*sin(j*t), 0, 0] g2 += [ 0, cos(i*x)*sin(j*t), 0] g3 += [ 0, 0, sin(i*x)*cos(j*t)] g = Matrix( [g1, g2, g3] ) B = Matrix(\ [[ D(x), 0, 0], [ 1/r*A, 0, 0], [ 1/r*D(t), 0, 0], [ 0, D(x), 0], [ 0, 1/r*A, 1/r*D(t)], [ 0, 1/r*D(t), D(x)-1/x], [ 0, 0, 1], [ 0, 1, 0]]) ans = evaluateExpr(B*g) 

Se ha creado una función print_to_file() para verificar rápidamente las expresiones grandes.

 import sympy import subprocess def print_to_file( guy, append=False ): flag = 'w' if append: flag = 'a' outfile = open(r'print.txt', flag) outfile.write('\n') outfile.write( sympy.pretty(guy, wrap_line=False) ) outfile.write('\n') outfile.close() subprocess.Popen( [r'notepad.exe', r'print.txt'] ) print_to_file( B*g ) print_to_file( ans, append=True ) 

Los operadores diferenciales no existen en el núcleo de SymPy, e incluso si existieran “multiplicación por un operador” en lugar de “aplicación de un operador” es un abuso de notación que no es compatible con SymPy.

[1] Otro problema es que las expresiones SymPy solo pueden comstackrse a partir de subclases de sympy.Basic , por lo que es probable que su class D simplemente sympy_expr+D(z) un error cuando se ingrese como sympy_expr+D(z) . Esta es la razón por la que (expression*D(z)) * (another_expr) falla. (expression*D(z)) no se puede construir.

Además, si el argumento de D no es un único Symbol no está claro qué espera de este operador.

Finalmente, diff(f(x), x) (donde f es una función desconocida simbólica) devuelve expresiones no evaluadas como usted observó simplemente porque cuando f es desconocido, no hay nada más que pueda ser devuelto sensiblemente. Más adelante, cuando sustituya expr.subs(f(x), sin(x)) se evaluará el derivado (en el peor de los casos, deberá llamar a expr.doit() ).

[2] No hay una solución elegante y breve para tu problema. Una forma que sugeriría para resolver su problema es anular el método __mul__ de Expr : en lugar de simplemente multiplicar los árboles de expresiones, comprobará si el árbol de expresiones de la izquierda contiene instancias de D y las aplicará. Obviamente, esto no se escala si desea agregar nuevos objetos. Este es un problema conocido desde hace mucho tiempo con el diseño de Sympy.

EDITAR: [1] es necesario simplemente para permitir la creación de expresiones que contengan D [2] es necesario para expresiones que contienen algo más que solo una D para que funcione.

Si desea que la multiplicación correcta funcione, deberá subclasificar solo el object . Eso hará que x*D caiga de nuevo a D.__rmul__ . No puedo imaginar que esto sea una prioridad alta, ya que los operadores nunca se aplican desde la derecha.

Hacer un operador que trabaje automáticamente siempre no es posible actualmente. Para trabajar realmente completamente, necesitarías http://code.google.com/p/sympy/issues/detail?id=1941 . Vea también https://github.com/sympy/sympy/wiki/Canonicalization (siéntase libre de editar esa página).

Sin embargo, podría crear una clase que funcione la mayor parte del tiempo utilizando las ideas de esa pregunta de stackoverflow, y para los casos que no maneja, escriba una función simple que pase por una expresión y aplique el operador donde no ha estado aplicado aún.

Por cierto, una cosa a considerar con un operador diferencial como “multiplicación” es que no es asociativa. A saber, (D*f)*g = g*Df , mientras que D*(f*g) = g*Df + f*Dg . Así que debes tener cuidado cuando hagas cosas que no “comen” una parte de una expresión y no todo. Por ejemplo, D*2*x daría 0 debido a esto. SymPy en todas partes asume que la multiplicación es asociativa, por lo que es probable que lo haga incorrectamente en algún momento.

Si eso se convierte en un problema, recomendaría descartar la aplicación automática, y simplemente trabajar con una función que se ejecute y aplique (que como he indicado anteriormente, la necesitará de todos modos).