差异运算符在Matrix模式中可用，在Python模块Sympy中-Java 学习之路

我们需要两个微分算子矩阵 [B] 和 [C] ，例如：

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]

这也可以应用于计算矢量场的卷曲，如：

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

要使用Sympy解决这个问题，必须创建以下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 )

当差分算子矩阵在左边相乘时，这个类就解决了 . 这里只有当要区分的函数已知时才执行 diff .

要在差分运算符矩阵乘以右侧时进行解决，必须按以下方式更改核心类 Expr 中的 __mul__ 方法：

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中应该有一个更好的原生解决方案来处理这个问题 . 有谁知道它可能是什么？

4 回答

此解决方案应用其他答案的提示和from here . D 运算符可以定义如下：

仅在从左侧乘以时才考虑，因此 D(t)*2*t**3 = 6*t**2 但 2*t**3*D(t) 什么都不做
D 使用的所有表达式和符号必须 is_commutative = False
使用 evaluateExpr() 在给定表达式的上下文中评估
从找到 D 操作符并将 mydiff() *应用于相应的右边部分的表达式从右向左移动

*：使用mydiff代替diff来允许创建更高阶D，如mydiff（D（t），t）= D（t，t）

D 中 __mul__() 内的 diff 仅供参考，因为在当前解决方案中， evaluateExpr() 实际上执行区分作业 . 创建了一个python mudule并保存为 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 中的顺序会影响结果，请参阅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近似 .

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() 函数以快速检查大表达式 .

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 )

回复于 2024-05-10T02:24:22+08:00

1

差异 operators 不存在于SymPy的核心中，即使它们存在"multiplication by an operator"而不是"application of an operator"，也是滥用SymPy不支持的符号 .

[1]另一个问题是SymPy表达式只能从 sympy.Basic 的子类构建，所以当输入为 sympy_expr+D(z) 时， class D 可能只会引发错误 . 这就是 (expression*D(z)) * (another_expr) 失败的原因 . (expression*D(z)) 无法构建 .

此外，如果 D 的参数不是单个 Symbol ，则不清楚您对此运算符的期望 .

最后， diff(f(x), x) （其中 f 是一个符号未知函数）返回一个未经评估的表达式，因为你观察的只是因为当 f 未知时，没有任何其他可以明智地返回 . 之后，当您替换 expr.subs(f(x), sin(x)) 时，将对导数进行求值（最坏的情况下，您可能需要调用 expr.doit() ） .

[2]对你的问题没有优雅的解决方案 . 我建议解决问题的方法是覆盖 Expr 的 __mul__ 方法：而不是仅仅将表达式树相乘，它将检查左表达式树是否包含 D 的实例，并且它将应用它们 . 显然，如果要添加新对象，则不会缩放 . 这是一个长期以来已知的问题设计问题 .

编辑：[1]只需要允许创建包含 D 的表达式 . [2]对于包含更多只有一个 D 工作的表达式是必要的 .

回复于 2024-05-10T02:24:22+08:00
1

如果你想要正确的乘法工作，你需要从 object 继承 . 这将导致 x*D 回落到 D.__rmul__ . 但是，我无法想象这是高优先级，因为操作员永远不会从右边应用 .

回复于 2024-05-10T02:24:22+08:00
3

目前无法实现自动运行的操作员 . 要真正完成工作，你需要http://code.google.com/p/sympy/issues/detail?id=1941 . 另请参阅https://github.com/sympy/sympy/wiki/Canonicalization（随意编辑该页面） .

但是，您可以使用来自stackoverflow问题的想法创建一个大部分时间都可以工作的类，对于它不处理的情况，编写一个通过表达式的简单函数，并将运算符应用于尚未处理的情况 . 应用了 .

顺便说一下，使用微分算子作为"multiplication"时要考虑的一点是它是非关联的 . 即， (D*f)*g = g*Df ，而 D*(f*g) = g*Df + f*Dg . 因此，当你做的事情不需要表达的某些部分而不是整个事情时，你需要小心 . 例如， D*2*x 会因此而给出 0 . SymPy无处不在假设乘法是关联的，因此在某些时候可能会错误地做到这一点 .

如果这成为一个问题，我建议转储自动应用程序，并只使用一个经过并应用它的功能（正如我上面提到的，你仍然需要） .

回复于 2024-05-10T02:24:22+08:00

差异运算符在Matrix模式中可用，在Python模块Sympy中

4 回答

相关问题