使用 Sympy 的方程系统的符号解，具有取决于符号的平凡解

Question

我有一个如下方程组：

对于这个特定的系统，我知道只有当 p1 == p2 时才存在非平凡的解决方案，即

.

但是，在一般情况下，如何使用 Sympy 确定这一点？

对于这个例子，我的实现如下：

from sympy import Matrix, symbols, pprint, lcm, latex
from sympy.solvers import solve_linear_system

top_matrix = Matrix.zeros(8,7)
p1 = symbols("p1")
p2 = symbols("p2")

top_matrix[0,0] = 1
top_matrix[0,1] = -1

top_matrix[1,1] = (1-p1)
top_matrix[1,2] = -1

top_matrix[2,2] = 1
top_matrix[2,4] = p2-1

top_matrix[3,1] = p1
top_matrix[3,3] = -1

top_matrix[4,3] = 1
top_matrix[4,4] = -p2

top_matrix[5,4] = 1
top_matrix[5,5] = -1

top_matrix[6,1] = -1
top_matrix[6,6] = 1

top_matrix[7,4] = -1
top_matrix[7,6] = 1

pprint(top_matrix)
vars = symbols("a1, a2, a3, a4, a5, a6, a7, a8")
print solve_linear_system(top_matrix, *vars)

结果是

None

如果我设置

p2 = p1

结果是

{a1: -1, a5: -1, a2: -1, a6: -1, a3: p1 - 1, a4: -p1}

有没有办法自动发现这个要求？

Answer 1

似乎 sympy 将 p1 和 p2 视为不相等，这意味着

top_matrix.nullspace()

[]

即矩阵系数top_matrix的齐次系统没有非平凡解。

这里有两个选择。第一种是将齐次系统视为具有未知变量的向量b和p2的元素，将p1视为固定参数。

b = sp.Matrix(sp.symbols('b1:8')) #import sympy as sp
sp.solve(top_matrix*b, (*b,p2))

[{b1: 0, b2: 0, b3: 0, b4: 0, b5: 0, b6: 0, b7: 0}, 
 {b1: b7, b2: b7, b3: b7*(-p1 + 1), b4: b7*p1, b5: b7, b6: b7, p2: p1}]

请注意，系统有两种解决方案。一个平凡的（任意p2），和一个非平凡的，它持有p2==p1（任意b7）。

第二个选项是如果系统A.T*A*b=0 有一个非平凡的解决方案，那么系统A*b=0 有一个非平凡的解决方案。如果A.T*A 的行列式为零，则后者是可能的。 A.T*A 的行列式等于

(top_matrix.T * top_matrix).det().factor()

6*(p1 - p2)**2

立即表明它必须持有 p1==p2 才能获得重要的解决方案。

Answer 2

在您的示例代码中，solve_linear_system 需要一个增强系统，即，如果右侧为零，则矩阵应声明为 Matrix.zeros(8,8)。通过此修改，您的代码将产生

{a3: 0, a1: 0, a5: 0, a7: 0, a6: 0, a2: 0, a4: 0}

这确实是一个解决方案，虽然不是很有趣...

要解决这个问题，可以明确要求将解决方案的一个组件标准化为 1。因此，如果您执行以下操作：

from sympy import Matrix, symbols, pprint, lcm, latex, solve

top_matrix = Matrix.zeros(8,7)
p1,p2 = symbols("p1, p2")

top_matrix[0,0] = 1
top_matrix[0,1] = -1

top_matrix[1,1] = (1-p1)
top_matrix[1,2] = -1

top_matrix[2,2] = 1
top_matrix[2,4] = p2-1

top_matrix[3,1] = p1
top_matrix[3,3] = -1

top_matrix[4,3] = 1
top_matrix[4,4] = -p2

top_matrix[5,4] = 1
top_matrix[5,5] = -1

top_matrix[6,1] = -1
top_matrix[6,6] = -1

top_matrix[7,4] = 1
top_matrix[7,6] = 1

pprint(top_matrix)

a1,a2,a3,a4,a5,a6,a7 = list(symbols("a1, a2, a3, a4, a5, a6, a7"))

B = Matrix([[1],[a2],[a3],[a4],[a5],[a6],[a7]])

C = top_matrix * B

print(solve(C, (a2,a3,a4,a5,a6,a7,p1,p2)))

并求解剩下的变量以及参数p1,p2，结果为：

[{a2: 1, a7: -1, a4: p2, a6: 1, a5: 1, p1: p2, a3: -p2 + 1}]

这确实是理想的解决方案。