我们需要微分运算符[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 中应该有一个更好的 native 解决方案来处理这个问题。 有人知道它可能是什么吗?
最佳答案
此解决方案应用了其他答案和 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 模块并保存为 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 )
关于python - 微分运算符可用于矩阵形式,在 Python 模块 Sympy 中,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/15463412/