python - Python-sympy 中的求解函数对于线性方程系统非常慢

Question

例如

from sympy import *

v11, v12, v21, v22 = symbols('v11, v12, v21, v22')
w11, w12, w21, w22 = symbols('w11, w12, w21, w22')
x11, x12, x21, x22 = symbols('x11, x12, x21, x22')
y11, y12, y21, y22 = symbols('y11, y12, y21, y22')
z11, z12, z21, z22 = symbols('z11, z12, z21, z22')

r11, r12, r21, r22 = symbols('r11, r12, r21, r22')

a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, \
a12, a13, a14, a15, a16, a17, a18 = symbols('a1:19')

# want to solve
b1, b2 = symbols('b1, b2')

# system of equations
eqn1  = a1*v11  + a2*v12  - b1       # a1 = 1
eqn2  = a1*v21  + a2*v22  - a3
eqn3  = a3*r11  + a4*r12  - a2
eqn4  = a3*r21  + a4*r22  - a5
eqn5  = a5*w11  + a6*w12  - a4   
eqn6  = a5*w21  + a6*w22  - a7
eqn7  = a7*r11  + a8*r12  - a6
eqn8  = a7*r21  + a8*r22  - a9
eqn9  = a9*x11  + a10*x12 - a8
eqn10 = a9*x21  + a10*x22 - a11
eqn11 = a11*r11 + a12*r12 - a10
eqn12 = a11*r21 + a12*r22 - a13
eqn13 = a13*y11 + a14*y12 - a12
eqn14 = a13*y21 + a14*y22 - a15
eqn15 = a15*r11 + a16*r12 - a14
eqn16 = a15*r21 + a16*r22 - a17
eqn17 = a17*z11 + a18*z12 - a16      # a18 = 0
eqn18 = a17*z21 + a18*z22 - b2

print('begin')
result = solve([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, \
                eqn8, eqn9, eqn10, eqn11, eqn12, eqn13, \
                eqn14, eqn15, eqn16, eqn17, eqn18], \
               [a2, a3, a4, a5, a6, a7, a8, a9, a10, \
                a11, a12, a13, a14, a15, a16, a17, b1, b2], \
                simplify=False, rational=False, manual=True)

b1_ans = result[b1]
b2_ans = result[b2]
print(b1_ans)

求解这些方程只需要几秒钟的 MATLAB，但使用 sympy 我无法得到结果。

有什么方法可以加快求解器的速度吗？如果不能，您能否推荐另一种使用 Python 求解方程组的方法？

提前致谢

Answer 1

由于您的符号比方程式多得多，因此您必须消除那些您不感兴趣的符号。大概您想要消除 want = set(var('a1:19')) 给定eqs = [Eq(a1, 1), Eq(a18, 0), eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn7, eqn8, eqn9, eqn10, eqn11, eqn12, eqn13, eqn14 的约束, eqn15, eqn16, eqn17, eqn18].

在递归地消除这些之后(求解只出现在 1 个方程中的变量然后从 want 中删除该变量，获得以下关系:

# recursively eliminate equations that have only one symbol of interest
sol = {}
while 1:
  was = len(sol)
  for ix, i in enumerate(eqs):
    if len(i.free_symbols & want) == 1:
        new = solve(i, *want, dict=True)
        assert len(new) == 1  # only doing this for single-solution equations
        sol.update(new[0])  # store in known values
  if len(sol) == was:
    break

方程式

0 a1 - 1
1 a18
2 a1*v11 + a2*v12 - b1
3 a1*v21 + a2*v22 - a3
4 -a2 + a3*r11 + a4*r12
5 a3*r21 + a4*r22 - a5
6 -a4 + a5*w11 + a6*w12
7 a5*w21 + a6*w22 - a7
8 -a6 + a7*r11 + a8*r12
9 a7*r21 + a8*r22 - a9
10 a10*x12 - a8 + a9*x11
11 a10*x22 - a11 + a9*x21
12 -a10 + a11*r11 + a12*r12
13 a11*r21 + a12*r22 - a13
14 -a12 + a13*y11 + a14*y12
15 a13*y21 + a14*y22 - a15
16 -a14 + a15*r11 + a16*r12
17 a15*r21 + a16*r22 - a17
18 -a16 + a17*z11 + a18*z12
19 a17*z21 + a18*z22 - b2

解决方案

{a1: 1, a18: 0, a2: (-a1*v11 + b1)/v12, a3: a1*v21 + a2*v22, a4: (a2 - a3*r11)/r12, a5: a3*r21 + a4*r22, a6: (a4 - a5*w11)/w12, a7: a5*w21 + a6*w22, a8: (a6 - a7*r11)/r12, a9: a7*r21 + a8*r22, a10: (a8 - a9*x11)/x12, a11: a10*x22 + a9*x21, a12: (a10 - a11*r11)/r12, a13: a11*r21 + a12*r22, a14: (a12 - a13*y11)/y12, a15: a13*y21 + a14*y22, a16: (a14 - a15*r11)/r12, a17: a15*r21 + a16*r22}

未解决/使用的方程式

{Eq(-a16 + a17*z11 + a18*z12, 0), Eq(a17*z21 + a18*z22 - b2, 0)}

可以为 b1 重新安排 a2 的解，并且可以为 b2 求解未解集中的第二个方程。

>>> solve(Eq(a17*z21 + a18*z22 - b2, 0), b2, dict=True)
[{b2: a17*z21 + a18*z22}]
>>> solve(s1[0][a2],b1, dict=True)
[{b1: a1*v11}]

剩下的方程表示一个必须满足的条件才能有解:Eq(-a16 + a17*z11 + a18*z12, 0)。

看看 MATLAB 的结果是什么会很有趣。