python - 如何在 SymPy 中扩展矩阵表达式？

Question

在 SymPy 中，我尝试执行矩阵乘法并在之后展开它。不过SymPy好像不支持矩阵表达式的展开。例如，这是矩阵的四阶龙格库塔 (RK4):

from sympy import init_session
init_session()
from sympy import *

A = MatrixSymbol('A', 3, 3)
x = MatrixSymbol('x', 3, 1)
dt = symbols('dt')

k1 = A*x
k2 = A*(x + S(1)/2*k1*dt)
k3 = A*(x + S(1)/2*k2*dt)
k4 = A*(x + k3*dt)
final = dt*S(1)/6*(k1 + 2*k2 + 2*k3 + k4)
final.expand()

产生结果

Traceback (most recent call last)
<ipython-input-38-b3ff67883c61> in <module>()
     12 final = dt*1/6*(k1+2*k2+2*k3+k4)
     13 
---> 14 final.expand()

AttributeError: 'MatMul' object has no attribute 'expand'

我希望表达式可以像标量变体一样展开:

A,x,dt = symbols('A x dt')
k1 = A*x
k2 = A*(x+k1*dt*S(1)/2)
k3 = A*(x+k2*dt*S(1)/2)
k4 = A*(x+k3*dt)
final = x+dt*(S(1)/6)*(k1+k2+k3+k4)
collect(expand((final)),x)

结果:

x*(A**4*dt**4/24 + A**3*dt**3/8 + A**2*dt**2/3 + 2*A*dt/3 + 1)

是否可以同样改变矩阵表达式？

nicoguaro 的回答消除了错误，但将整个表达式扩展为一个矩阵。正如标量示例所示，这不是我要找的。

Answer 1

Matrix(final) 使用您的个人方程式创建和显式矩阵。将它们留在矩阵中可能很方便，这样您对一个条目所做的任何事情都可以对所有条目进行。为此，将您想要执行的操作定义为 applyfunc 的函数参数:

>>> ex = Matrix(final)
>>> ex = ex.applyfunc(expand)
>>> ex = ex.applyfunc(lambda i: collect(i, dt))
...

为了打印经济，我使用具有紧凑符号条目的矩阵来计算相同的值，然后在简化矩阵上运行 cse 以给出:

>>> A = Matrix(3, 3, var('a:3:3'))
>>> x = Matrix(3, 1, var('x:3'))
>>> dt = symbols('dt')
>>> k1 = A*x
>>> k2 = A*(x + S(1)/2*k1*dt)
>>> k3 = A*(x + S(1)/2*k2*dt)
>>> k4 = A*(x + k3*dt)
>>> final = dt*S(1)/6*(k1 + 2*k2 + 2*k3 + k4)
>>> eqi = Matrix(final)
>>> print(cse(eqi.applyfunc(simplify)))
([(x3, 4*x0), (x4, a01*x1), (x5, a02*x2), (x6, dt*(a00*x0 + x4 + x5) +
2*x0), (x7, a00*x6), (x8, a11*x1), (x9, a12*x2), (x10, dt*(a10*x0 + x8
+ x9) + 2*x1), (x11, a01*x10), (x12, a21*x1), (x13, a22*x2), (x14,
dt*(a20*x0 + x12 + x13) + 2*x2), (x15, a02*x14), (x16, dt*(x11 + x15 +
x7) + x3), (x17, a00*x16), (x18, 4*x1), (x19, a10*x6), (x20, a11*x10),
(x21, a12*x14), (x22, dt*(x19 + x20 + x21) + x18), (x23, a01*x22),
(x24, 4*x2), (x25, a20*x6), (x26, a21*x10), (x27, a22*x14), (x28,
dt*(x25 + x26 + x27) + x24), (x29, a02*x28), (x30, dt*(x17 + x23 +
x29) + x3), (x31, a10*x16), (x32, a11*x22), (x33, a12*x28), (x34,
dt*(x31 + x32 + x33) + x18), (x35, a20*x16), (x36, a21*x22), (x37,
a22*x28), (x38, dt*(x35 + x36 + x37) + x24), (x39, dt/24)], [Matrix([
[   x39*(a00*x3 + a00*x30 + a01*x34 + a02*x38 + 4*x11 + 4*x15 + 2*x17
+ 2*x23 + 2*x29 + 4*x4 + 4*x5 + 4*x7)], [  x39*(a10*x3 + a10*x30 +
a11*x34 + a12*x38 + 4*x19 + 4*x20 + 4*x21 + 2*x31 + 2*x32 + 2*x33 +
4*x8 + 4*x9)], [x39*(a20*x3 + a20*x30 + a21*x34 + a22*x38 + 4*x12 +
4*x13 + 4*x25 + 4*x26 + 4*x27 + 2*x35 + 2*x36 + 2*x37)]])])

非交换表达式的有限支持也可用，在这种情况下可以提供帮助:

>>> A, x = symbols("A x", commutative=False)
>>> dt = symbols('dt')
>>> k1 = A*x
>>> k2 = A*(x + S(1)/2*k1*dt)
>>> k3 = A*(x + S(1)/2*k2*dt)
>>> k4 = A*(x + k3*dt)
>>> final = dt*S(1)/6*(k1 + 2*k2 + 2*k3 + k4)
>>> final.expand()
dt**4*A**4*x/24 + dt**3*A**3*x/6 + dt**2*A**2*x/2 + dt*A*x
>>> factor(_)
dt*A*(dt**3*A**3/24 + dt**2*A**2/6 + dt*A/2 + 1)*x

但并非所有简化例程都是 nc 感知的(这是一个已知问题):

>>> collect(final,x)
Traceback (most recent call last):
...
AttributeError: Can not collect noncommutative symbol