可在Python模块Sympy中以矩阵形式使用的微分运算符 [英] Differential Operator usable in Matrix form, in Python module Sympy

问题描述

我们需要两个差分运算符[B]和[C]的矩阵，例如:

We need two matrices of differential operators [B] and [C] such as:

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]

这也可以用于计算矢量场的卷曲度，例如:

This could also be applied to calculate the curl of a vector field like:

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

要使用Sympy解决此问题，必须创建以下Python类:

To solve this using Sympy the following Python class had to be created:

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 )

仅此类可解决差分运算符的矩阵在左侧相乘的情况.这里diff仅在已知要区分的函数时执行.

This class alone solves when the matrix of differential operators is multiplying on the left. Here diff is executed only when the function to be differentiated is known.

要解决变数运算符矩阵在右侧相乘的问题，必须以下列方式更改核心类Expr中的__mul__方法:

To workaround when the matrix of differential operators is multiplying on the right, the __mul__ method in the core class Expr had to be changed in the following way:

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)
    #...

它工作得很好，但是Sympy中应该有一个更好的本机解决方案来处理此问题. 有人知道这可能是什么吗?

It works pretty well, but there should be a better native solution in Sympy to handle this. Does anybody know what it might be?

推荐答案

此解决方案适用于其他答案和

This solution applies the tips from the other answers and from here. The D operator can be defined as follows:

• 仅当从左侧乘以时考虑，因此D(t)*2*t**3 = 6*t**2但2*t**3*D(t)无效
• 与D一起使用的所有表达式和符号都必须具有is_commutative = False
使用evaluateExpr()在给定表达式的上下文中评估
• • 沿着表达式从右到左查找D运算符，并将mydiff() *应用于相应的右侧部分
  • considered only when multiplied from the left, so that D(t)*2*t**3 = 6*t**2 but 2*t**3*D(t) does nothing
  • all the expressions and symbols used with D must have is_commutative = False
  • is evaluated in the context of a given expression using evaluateExpr()
    • which goes from the right to the left along the expression finding the D opperators and applying mydiff()* to the corresponding right portion

    *: _{mydiff代替diff用于允许创建更高阶的D，例如mydiff(D(t), t) = D(t,t)}

    *:_{mydiff is used instead of diff to allow a higher order D to be created, like mydiff(D(t), t) = D(t,t)}

    D中__mul__()内部的diff仅作参考，因为在当前解决方案中，evaluateExpr()实际上完成了微分工作.创建了一个python模版并将其另存为d.py.

    The diff inside __mul__() in D was kept for reference only, since in the current solution the evaluateExpr() actually does the differentiation job. A python mudule was created and saved as 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)
    

    示例1:向量场的卷曲.请注意，使用commutative=False定义变量很重要，因为它们在Mul().args中的顺序会影响结果，请参见

    Example 1: curl of a vector field. Note that it is important to define the variables with commutative=False since their order in Mul().args will affect the results, see this other question.

    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]
    

    示例2:结构分析中使用的典型Ritz近似值.

    Example 2: Typical Ritz approximation used in structural analysis.

    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)
    

    已创建print_to_file()函数以快速检查大表达式.

    A print_to_file() function has been created to quickly check big expressions.

    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 )
    

    这篇关于可在Python模块Sympy中以矩阵形式使用的微分运算符的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！