如何使用sympy简化矩阵中的分数？

Question

我正在使用sympy来找到一个矩阵的逆。我遇到了一个问题。当我计算矩阵A的逆时，我得到了一个带有分数的矩阵；我的意思是

>> import sympy
>> from sympy import pprint
>> from sympy.abc import *
>> import sys
>> sys.displayhook = pprint
>> from sympy.matrices import *
>> A = Matrix([[a, b],[c, d]])
>> B = A.inv()
>> B
>> [1       b*c           -b     ]
>> [- + ------------  -----------]
>> [a    2 /    b*c\    /    b*c\]
>> [    a *|d - ---|  a*|d - ---|]
>> [       \     a /    \     a /]
>> [                             ]
>> [      -c               1     ]
>> [  -----------       -------  ]
>> [    /    b*c\           b*c  ]
>> [  a*|d - ---|       d - ---  ]
>> [    \     a /            a   ]
>> B*A
>> [  /1       b*c     \       b*c        /1       b*c     \       b*d    ]
>> [a*|- + ------------| - -----------  b*|- + ------------| - -----------]
>> [  |a    2 /    b*c\|     /    b*c\    |a    2 /    b*c\|     /    b*c\]
>> [  |    a *|d - ---||   a*|d - ---|    |    a *|d - ---||   a*|d - ---|]
>> [  \       \     a //     \     a /    \       \     a //     \     a /]
>> [                                                                      ]
>> [                                             d          b*c           ]
>> [                0                         ------- - -----------       ]
>> [                                              b*c     /    b*c\       ]
>> [                                          d - ---   a*|d - ---|       ]
>> [                                               a      \     a /       ]

我想获取下一个矩阵。

>> I = Matrix([
>> [1, 0],
>> [0, 1]])

我的问题是矩阵 A*B 或 B*A。我想简化矩阵 A*B 并得到 I。我尝试使用 simplify() 但没有成功。

- Yakz

---

简化是什么意思？这是一些技术术语还是什么？我的第二个问题是，你要证明什么？ - marmeladze

好的，B是A的逆矩阵，它们的乘积应该是矩阵单位I。在代码的乘积中，我得到了一个巨大的矩阵AB（或BA），但我希望这个矩阵的条目为零（0）。sympy不会简化条目（1,2），（1,1）和（2,1）。您可以查看此乘积的结果。请原谅我的英语，我不太擅长。感谢您的时间。 - Yakz

我得到了与你相同的AB结果，但是simplify(AB)会给出单位矩阵。我在Windows 7上使用最新的Anaconda 64笔记本。 - brian

Answer 1

您可以使用applyfunc将simplify函数应用于矩阵的每个单元格，如下所示：

>>> (B*A).applyfunc(simplify)
[1  0]
[    ]
[0  1]

- Adrien

---

我正在尝试这样做，但是我得到了“名称'simplify'未定义”的错误。 - Yakz

是的，但没有什么。我编写了一个按条目执行此操作的函数。 - Yakz

Answer 2

先不要考虑Python和Sympy。用纸和笔专注于找到矩阵的逆。

对于一个A = [[a, b], [c,d]]矩阵，我们计算逆矩阵A^-1如下：

(1/D)*[[d, -b],[-c, a]]。这里D是矩阵A的行列式（1/ad-bc）

这个（A^-1）等于[[d/D, -b/D][-c/D, a/D]]

让我们从第一行的第一个元素开始，并按照我所做的操作进行。对我来说，它们实际上没有任何意义，但这是Sympy的工作方式 :) 然后将此过程应用于其他矩阵元素。

=> d/D 
d/(a*d-b*c) 
a*d/(d*a^2 - a*b*c)
(a*d-b*c+b*c)/a^2*(d-b*c/a) 
(a*d - a*b*c/a + b*c)/a^2*(d-b*c/a)
(a*(d-b*c/a) + b*c)/a^2*(d-b*c/a)
a*(d-b*c/a)/a^2*(d-b*c/a) + b*c/a^2*(d-b*c/a)
1/a + b*c/a^2*(d-b*c/a) [this is how sympy outputs]


>>> A = Matrix([[a,b],[c,d]])
>>> B = A**-1 #same as B = A.inv()
>>> B[0]
1/a + b*c/(a**2*(d - b*c/a))

现在，让我们来看一下sympy A*B的输出是什么。

>>> N = A*B
>>> N
Matrix([
[a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a)),                                   0],
[c*(1/a + b*c/(a**2*(d - b*c/a))) - c*d/(a*(d - b*c/a)), d/(d - b*c/a) - b*c/(a*(d - b*c/a))]])
>>> pprint(N)
⎡  ⎛1       b⋅c     ⎞       b⋅c                           ⎤
⎢a⋅⎜─ + ────────────⎟ - ───────────            0          ⎥
⎢  ⎜a    2 ⎛    b⋅c⎞⎟     ⎛    b⋅c⎞                       ⎥
⎢  ⎜    a ⋅⎜d - ───⎟⎟   a⋅⎜d - ───⎟                       ⎥
⎢  ⎝       ⎝     a ⎠⎠     ⎝     a ⎠                       ⎥
⎢                                                         ⎥
⎢  ⎛1       b⋅c     ⎞       c⋅d         d          b⋅c    ⎥
⎢c⋅⎜─ + ────────────⎟ - ───────────  ─────── - ───────────⎥
⎢  ⎜a    2 ⎛    b⋅c⎞⎟     ⎛    b⋅c⎞      b⋅c     ⎛    b⋅c⎞⎥
⎢  ⎜    a ⋅⎜d - ───⎟⎟   a⋅⎜d - ───⎟  d - ───   a⋅⎜d - ───⎟⎥
⎣  ⎝       ⎝     a ⎠⎠     ⎝     a ⎠       a      ⎝     a ⎠⎦

它不会直接评估为eye(2)，但如果您取元素并简化它们，您会发现这个混乱的矩阵实际上是一个2x2的单位矩阵。

检查它的Pythonic方法（已知）：

>>> N[0]
a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a))
>>> N[1]
0
>>> N[3]
d/(d - b*c/a) - b*c/(a*(d - b*c/a))
>>> N[2]
c*(1/a + b*c/(a**2*(d - b*c/a))) - c*d/(a*(d - b*c/a))

>>> def will_evaluate_one(a,b,c,d):
...    return a*(1/a + b*c/(a**2*(d - b*c/a))) - b*c/(a*(d - b*c/a)) #N[0]     
...
>>> will_evaluate_one(1,2,3,9)
1
>>> will_evaluate_one(1,2,3,19)
1
>>> will_evaluate_one(1,2,23,19)
1
>>> will_evaluate_one(1,12,23,19)
1

>>> def will_also_evaluate_one(a,b,c,d):
...     return d/(d - b*c/a) - b*c/(a*(d - b*c/a)) #N[1] 
... 
>>> will_also_evaluate_one(2,4,5,6)
1
>>> will_also_evaluate_one(2,4,15,6)
1
>>> will_also_evaluate_one(2,14,15,6)
1
>>> will_also_evaluate_one(12,14,15,6)
1

注意：我刚刚意识到，sympy使用解析反演公式。请参见此处：https://en.wikipedia.org/wiki/Helmert%E2%80%93Wolf_blocking。

- marmeladze