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

Question

I have an equation system like the following:

For this specific system, I know that a nontrivial solution(s) only exists if p1 == p2, which is

.

However, how can I determine this in the general case using Sympy?

For this example, my implementation is as follows:

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)

The result is

None

If I set

p2 = p1

the result is

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

Is there some way to discover this requirement automatically?

Answer 1

It appears that sympy treats p1 and p2 as non-equal, which means that

top_matrix.nullspace()

[]

that is, the homogeous system with matrix coefficients top_matrix has no non-trivial solution.

You have two options, here. The first is to treat the homogeneous system as having unknown variables both the elements of vector b and p2, treating p1 as a fixed parameter.

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}]

Note that the system has two solutions. The trivial one (with arbitrary p2), and a non-trivial one, where it holds p2==p1 (with arbitrary b7).

The second option is to realize that the system A*b=0 has a non-trivial solution if the system A.T*A*b=0 has a non-trivial solution. The latter is possible if the determinant of A.T*A is zero. The determinant of A.T*A equals

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

6*(p1 - p2)**2

immediately revealing that it must hold p1==p2 in order to have a non-trivial solution.

Answer 2

In your example code, solve_linear_system expects an augmented system, i.e., if the right hand side is zero, the matrix should be declared as Matrix.zeros(8,8). With this modification, your code yields

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

which is indeed a solution, although not the interesting one...

To remedy this, one can explicitly request that one component of the solution should be normalized to, say, 1. So if you do something like the following:

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)))

and solve for the remaining variables as well as the parameters p1,p2, the result is:

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

which is indeed the desired solution.