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

Symbolic solution of equation system using Sympy with trivial solutions depending on symbols

我有一个方程式如下:

对于这个特定的系统,我知道只有当 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}

有什么方法可以自动发现这个要求吗?

似乎 sympy 将 p1p2 视为不相等,这意味着

top_matrix.nullspace()

[]

即矩阵系数为top_matrix的齐次系统无非平凡解

这里有两个选择。第一种是将齐次系统视为具有向量bp2两个元素的未知变量,将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 才能得到一个非平凡的解决方案。

在您的示例代码中,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}]

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