# Sage wont solve simultaneous equation

I tried asking this question earlier, but my profile still says I have asked 0 questions. I have tried solving this problem in different ways, but this is the closest I have gotten;

var('a b c d e f g h i')
X=matrix(3,3,[[0,1,0],[0,0,1],[1,0,0]])
P=matrix(3,3,[[a,b,c],[d,e,f],[g,h,i]])
Pdagger=P.transpose()
Xdagger=X.transpose()
Q=X*Pdagger
B=Q.solve_left(X)
eqns=[]
xlist=[]
for i in range(0,3):
for j in range(0,3):
eqns.append(B[i][j])
xlist.append(X[i][j])
final=[]
for k in range(0,9):
final.append(eqns[k]==xlist[k])
show(final)
sol=solve(eqns,[a,b,c,d,e,f,g,h,i])
show(sol)


I know there is a solution to this problem, and it is the identity, $PXP^{\dagger}=X$.

But the output:

Error in lines 19-19
Traceback (most recent call last):
File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1013, in execute
exec compile(block+'\n', '', 'single') in namespace, locals
File "", line 1, in <module>
File "/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/symbolic/relation.py", line 976, in solve
raise TypeError("%s is not a valid variable." % repr(i))
TypeError: 2 is not a valid variable.


There is also probably a much more efficient way of specifying the problem, but I don't understand why this wont work. I would like to extend this problem to larger matrices as well, so autonomy is welcomed. However, I would appreciate any advice just on getting this error fixed.

edit retag close merge delete

This seems to trigger a bug. I'll need some time to follow it (not so soon).

( 2018-04-26 12:47:39 -0500 )edit

Sort by ยป oldest newest most voted

It is hard to understand what happens behind, but the following works:

vars = var('a b c d e f g h u')

X = matrix(3, 3, [[0, 1, 0], [0, 0, 1], [1, 0, 0]])
P = matrix(3, 3, [[a, b, c], [d, e, f], [g, h, u]])

Pdagger = P.transpose()
Xdagger = X.transpose()    # not needed

Q = X * Pdagger
B = Q.solve_left(X)

R     = range(0,3)
eqns  = [ B[i][j] for i in R for j in R ]
xlist = [ X[i][j] for i in R for j in R ]
final = [ B[i][j] == X[i][j] for i in R for j in R ]

sols = solve(final, vars)
# show(sol)
for sol in sols:
print sol


• when programming, please use a decent spacing, else nobody will see an error when an obvious one exists. "pep8".
• if i is a variable... ok, no problem we overwrite that $\sqrt {-1}$, then please do not redefine it in a loop as a number. (This would have been a second alarm while reusing i...) Nobody does this, so this atypical error cannot be seen at a first glance. In our case this brings the error as follows. After the last i-loop this "global variable" has the value $2$, and we try to solve a system of equations, where one of the equation variables is the constant $2$. The error message is also misleading.
• give a minimal code showing the error, only needed lines of code.
• use list comprehension.
• do not call a list eqns, when it does not collect equations, but rather entries of a matrix. The solve of these "equations" will be hardly detected with bare eyes as an error.
• try to give pointed prints. The show command delivers a mess, is not a good control of the operations.
• the two definitions of eqns and xlist are not needed. The wording for xlist is also not proper. Something like X_entries would reflect the provenience, but we do not need this new variable, the matrix X already has the information.

Please understand that these "critics" are just kind advices, at any rate not offensive, in my case, knowing these as i started programming in python would have been a real benefit.

more

Let's be lazy (and avoid identifier clashes)

sage: V=var("v", n=9)
sage: P=matrix(V,nrows=3)
sage: X=matrix(SR,3,3,[[0,1,0],[0,0,1],[1,0,0]])
sage: Pdagger=P.transpose()
sage: Xdagger=X.transpose()
sage: Q=X*Pdagger
sage: B=Q.solve_left(X)


Let's be lazy again : don't build a list of equations representing B==X, but a list representing the elements of B-X

sage: E=[B[r][c]-X[r][c] for c in range(3) for r in range(3)]


Look (lazily) for solutions

sage: S=solve(E,V)


A loook at the solutions show that they are rational expressions. Are they acceptable ? This entails that their denominators must not be null. But :

sage: [[e.denominator().subs(s) for s in S] for e in E]
[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1]]


Check the only acceptable one :

sage: [e.subs(S[16]) for e in E]
[0, 0, 0, 0, 0, 0, 0, 0, 0]


And, indeed :

sage: matrix([v.subs(S[16]) for v in V],nrows=3)==X
True


Note that Sympy does not seem to be able to solve this system :

sage: sympy.solve(map(lambda e:e._sympy_(), E),map(lambda v:v._sympy_(), V))
[]


On the other hand, Mathematica is able to solve the system and does not gives "spurious" solutions. But expressing the exuations is harder in Mathematica than in Sage : that's why I keep the first part in Sage :

sage: V=var("v",n=9)
sage: P=matrix(SR,V,nrows=3);P
[v0 v1 v2]
[v3 v4 v5]
[v6 v7 v8]
sage: X=matrix(SR,[[0,1,0],[0,0,1],[1,0,0]]);X
[0 1 0]
[0 0 1]
[1 0 0]
sage: Q=X*(P.transpose())
sage: B=Q.solve_left(X)
sage: E=[B[l][c]==X[l][c] for c in range(3) for l in range(3)];E
[(v2*v7/v1 - v8)*(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) + v2/(v1*(v2*v4/v1 - v5)) == 0,
-(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) == 0,
-(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))*(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) + 1/v1 - v2*v4/(v1^2*(v2*v4/v1 - v5)) == 1,
-1/(v2*v4/v1 - v5) - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/((v2*v4/v1 - v5)^2*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 1,
(v3 - v0*v4/v1)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 0,
(v3 - v0*v4/v1)*(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) + v4/(v1*(v2*v4/v1 - v5)) == 0,
-(v2*v7/v1 - v8)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 0,
(1/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 1,
(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) == 0]
sage: mathematica.Solve(mathematica.And(*E), V)
{{v0 -> 0, v1 -> 1, v2 -> 0, v3 -> 0, v4 -> 0, v5 -> 1, v6 -> 1, v7 -> 0,
v8 -> 0}}

more