用符号替换子表达式？
我有一个3x3的矩阵，我想计算它的逆矩阵。为了让逆矩阵看起来更清晰，我需要把一些重复出现的部分用新的符号替换掉。这样做是因为这些部分出现了很多次。请问，我能不能让sympy尽量找到这些重复的部分并替换掉呢？我试过以下的方法，但没有成功：
from sympy import *
Ex, Ez, nuxy, nuxz = symbols('E_x E_z nu_xy nu_xz')
# compliance matrix for cross-anisotropic material
compl = Matrix([[1/Ex, -nuxy/Ex, -nuxz/Ez],
[-nuxy/Ex, 1/Ex, -nuxz/Ez],
[-nuxz/Ex, -nuxz/Ex, 1/Ez]])
# stiffness matrix
stiff = compl.inv()
# symbols I want to introduce
m, e = symbols('m e')
meSubs = {Ex/Ez: e, (1 - nuxy - 2*e*nuxz**2): m} # instead of these subexpressions
# stiff.simplify() returns None, is that a bug? that's why I apply simplify together with subs here:
stiff.applyfunc(lambda x: simplify(x.subs(meSubs)))
print stiff
我使用的是sympy 0.6.7（如果需要的话，我可以升级）。
编辑：
我把版本升级到了0.7.1-git（具体是cf9c01f8f9b4b749a7f59891f546646e4b38e580），并运行了（感谢@PreludeAndFugue的建议）：
from sympy import *
Ex,Ez,nuxy,nuxz,m=symbols('E_x E_z nu_xy nu_xz m')
compl=Matrix([[1/Ex,-nuxy/Ex,-nuxz/Ez],[-nuxy/Ex,1/Ex,-nuxz/Ez],[-nuxz/Ex,-nuxz/Ex,1/Ez]])
stiff=compl.inv()
stiff.simplify()
stiff.subs({-nuxy-2*nuxz**2+1:m}) # tried other rearrangements of the equality, as well, same result.
stiff.applyfunc(lambda x: together(expand(x)))
pprint(stiff)
得到了
⎡ ⎛ 2 ⎞ ⎛ 2⎞ ⎤
⎢ Eₓ⋅⎝ν_xz - 1⎠ -Eₓ⋅⎝-ν_xy - ν_xz ⎠ Eₓ⋅ν_xz ⎥
⎢ ────────────────────────────────── ──────────────────────────────────── ───────────────────⎥
⎢ 2 2 2 2 2 2 2 ⎥
⎢ ν_xy + 2⋅ν_xy⋅ν_xz + 2⋅ν_xz - 1 - ν_xy - 2⋅ν_xy⋅ν_xz - 2⋅ν_xz + 1 -ν_xy - 2⋅ν_xz + 1⎥
⎢ ⎥
⎢ ⎛ 2⎞ ⎛ 2 ⎞ ⎥
⎢ -Eₓ⋅⎝-ν_xy - ν_xz ⎠ Eₓ⋅⎝ν_xz - 1⎠ Eₓ⋅ν_xz ⎥
⎢──────────────────────────────────── ────────────────────────────────── ───────────────────⎥
⎢ 2 2 2 2 2 2 2 ⎥
⎢- ν_xy - 2⋅ν_xy⋅ν_xz - 2⋅ν_xz + 1 ν_xy + 2⋅ν_xy⋅ν_xz + 2⋅ν_xz - 1 -ν_xy - 2⋅ν_xz + 1⎥
⎢ ⎥
⎢ E_z⋅ν_xz E_z⋅ν_xz E_z⋅(ν_xy - 1) ⎥
⎢ ─────────────────── ─────────────────── ────────────────── ⎥
⎢ 2 2 2 ⎥
⎣ -ν_xy - 2⋅ν_xz + 1 -ν_xy - 2⋅ν_xz + 1 ν_xy + 2⋅ν_xz - 1 ⎦
嗯，那为什么"-ν_xy - 2⋅ν_xz² + 1"没有被替换成m呢？
3 个回答
0
确实会被替换，但是 subs 在矩阵上并不能直接修改。很遗憾，applyfunc 也不能直接修改。我得到的结果是
In [10]: pprint(stiff.subs({-nuxy-2*nuxz**2+1:m}))
⎡ ⎛ 2 ⎞ ⎛ 2⎞ ⎤
⎢ Eₓ⋅⎝nu_xz - 1⎠ -Eₓ⋅⎝-nu_xy - nu_xz ⎠ Eₓ⋅nu_xz ⎥
⎢ ────────────────────────────────────── ──────────────────────────────────────── ──────── ⎥
⎢ 2 2 2 2 2 2 m ⎥
⎢ nu_xy + 2⋅nu_xy⋅nu_xz + 2⋅nu_xz - 1 - nu_xy - 2⋅nu_xy⋅nu_xz - 2⋅nu_xz + 1 ⎥
⎢ ⎥
⎢ ⎛ 2⎞ ⎛ 2 ⎞ ⎥
⎢ -Eₓ⋅⎝-nu_xy - nu_xz ⎠ Eₓ⋅⎝nu_xz - 1⎠ Eₓ⋅nu_xz ⎥
⎢──────────────────────────────────────── ────────────────────────────────────── ──────── ⎥
⎢ 2 2 2 2 2 2 m ⎥
⎢- nu_xy - 2⋅nu_xy⋅nu_xz - 2⋅nu_xz + 1 nu_xy + 2⋅nu_xy⋅nu_xz + 2⋅nu_xz - 1 ⎥
⎢ ⎥
⎢ E_z⋅nu_xz E_z⋅nu_xz -E_z⋅(nu_xy - 1)⎥
⎢ ───────── ───────── ────────────────⎥
⎣ m m m ⎦
目前有计划让矩阵默认是不可变的，然后让可变矩阵在所有操作中都能完全就地修改。具体可以查看这个链接：https://code.google.com/p/sympy/issues/detail?id=3410。不过这项计划还没有实现。
0
我不太确定使用0.6.7版本是否有问题，但建议更新到0.7.1版本。
当我查看stiff时，发现meSubs中的替换似乎没什么用。在创建stiff之后，我做了以下操作：
stiff.simplify()
stiff = stiff.subs({2*nuxz**2: 1 - nuxy - m})
stiff = stiff.applyfunc(lambda x: together(expand(x)))
pprint(stiff)
输出结果还不错：
[ / 2 \ / 2\ ]
[ E_x*\nu_xz - 1/ E_x*\nu_xy + nu_xz / E_x*nu_xz ]
[ ---------------- -------------------- --------- ]
[ m*(-nu_xy - 1) m*(nu_xy + 1) m ]
[ ]
[ / 2\ / 2 \ ]
[E_x*\nu_xy + nu_xz / E_x*\nu_xz - 1/ E_x*nu_xz ]
[-------------------- ---------------- --------- ]
[ m*(nu_xy + 1) m*(-nu_xy - 1) m ]
[ ]
[ E_z*nu_xz E_z*nu_xz E_z*(-nu_xy + 1)]
[ --------- --------- ----------------]
[ m m m ]
展开：http://docs.sympy.org/0.7.1/modules/core.html?highlight=expand#sympy.core.function.expand
1
你们的替换都失败了，因为你们想要的模式在矩阵中不存在。我们一个一个来看：
- 没有Ex/Ez的比率；
x.subs(x/y, z)没有变化。 - 在
1 - nuxy - 2*e*nuxz**2中有一个多余的e，所以这个表达式也无法匹配。
@asmeurer展示了字面上的1 - nuxy - 2*nuxz**2可以被替换，但它也出现在一些分母中。更复杂的替换可以通过检查模式是否能整除一个表达式来完成。
让我们来写一个函数，来进行这个替换：
>>> from sympy import *
>>> t = 1 - nuxy - 2*e*nuxz**2
>>> def do(x):
... w, r = div(x, t)
... if not r:
... return m*w
... return x
现在我们将把这个应用到矩阵的每个元素上：
>>> stiff.applyfunc(lambda x: factor_terms(bottom_up(x, do)))
Matrix([
[ -E_x*(nu_xz**2 - 1)/(m*(nu_xy + 1)), E_x*(nu_xy + nu_xz**2)/(m*(nu_xy + 1)), E_x*nu_xz/m],
[E_x*(nu_xy + nu_xz**2)/(m*(nu_xy + 1)), -E_x*(nu_xz**2 - 1)/(m*(nu_xy + 1)), E_x*nu_xz/m],
[ E_z*nu_xz/m, E_z*nu_xz/m, -E_z*(nu_xy - 1)/m]])
对于复杂的表达式或矩阵，有时候使用cse可以帮助我们更好地了解结构：
>>> cse(stiff)
([(x0, nu_xz**2), (x1, E_x*x0), (x2, 2*x0), (x3, x2 - 1), (x4, 1/(nu_xy**2 + nu_xy*x2 + x3)), (x5, x4*(-E_x + x1)), (x6, x4*(-E_x*nu_xy - x1)), (x7, 1/(nu_xy + x3)), (x8, nu_xz*x7), (x9, -E_x*x8), (x10, -E_z*x8)], [Matrix([
[ x5, x6, x9],
[ x6, x5, x9],
[x10, x10, x7*(E_z*nu_xy - E_z)]])])