# How to stop Sage from finding erroneous complex roots?

Consider the matrix,

A = matrix ( [0,1,w^a,1],[1,0,1,w^(k-b)],[w^(k-a),1,0,w^(k-c)],[1,w^b,w^c,0])


where w is the m^th of the k^th roots of unity, w = exp((2*pi*I*m )/k for some k a positive integer, and 1 <= m <= (k-1) and 1<= a,b,c <= (k-1)

Then the characteristic polynomial of the above matrix is, p(x) = x^4 - 6*x^2 -x *(w^(a-c) + w^(c-a) + w^b + w^(-b) + w^(b-c) + w^(c-b) + w^a + w^(-a)) + (3 -w^c - w^(-c) - w^(a+b-c) - w^(-a-b+c) - w^(a-b) - w^(-a+b) )

• Is there a way to get sage to be able to calculate the above characteristic polynomial?

Now I try getting roots of the above by doing,

g(x)=real_part(p(x)).simplify()
g.solve(x)

• Now Sage seems to be generically detecting complex eigenvalues as roots of g! (I tried on say k=6, m=a=b=c=1) This can't happen since its a characteristic polynomial of a Hermitian matrix!

How to get across this trouble?

edit retag close merge delete

Sort by » oldest newest most voted

The main issue is that you deal with symbolic expressions (where computations are hard, for example checking that some symbolic expression is equal to zero is undecidable), not algebraic or numerical ones:

sage: k=6 ; m=a=b=c=1
sage: w = exp((2*pi*I*m )/k)
sage: A = matrix([[0,1,w^a,1],[1,0,1,w^(k-b)],[w^(k-a),1,0,w^(k-c)],[1,w^b,w^c,0]])
sage: A
[           0            1 e^(1/3*I*pi)            1]
[           1            0            1 e^(5/3*I*pi)]
[e^(5/3*I*pi)            1            0 e^(5/3*I*pi)]
[           1 e^(1/3*I*pi) e^(1/3*I*pi)            0]
sage: A.parent()
Full MatrixSpace of 4 by 4 dense matrices over Symbolic Ring


You can compute the characteristic polynomial:

sage: P = A.characteristic_polynomial()
sage: P
x^4 - 6*x^2 - 6*x - 1


But it has symbolic coefficients, so it is not able to find its roots:

sage: P.parent()
Univariate Polynomial Ring in x over Symbolic Ring
sage: P.roots()
NotImplementedError:


What you should do is to specify a better ring for the coefficients. You will see that the roots depend on the ring:

There is one integer root with multiplicity 1:

sage: Q = P.change_ring(ZZ)
sage: Q.roots()
[(-1, 1)]


There is one rational root with multiplicity 1:

sage: Q = P.change_ring(QQ)
sage: Q.roots()
[(-1, 1)]


They are 4 algebraic roots, each with multipicity 1.

sage: Q = P.change_ring(QQbar)
sage: Q.roots()
[(-1.655442381549831?, 1),
(-1, 1),
(-0.2107558809591918?, 1),
(2.866198262509023?, 1)]


You can also work in inexact rings, to get floating point approximations, for example in the Real Double Field:

sage: Q = P.change_ring(RDF)
sage: Q.roots()
[(-1.6554423815498323, 1),
(-0.9999999999999978, 1),
(-0.21075588095919173, 1),
(2.866198262509024, 1)]


Note that there is a huge difference between the last two blocks, the one with QQbar gives you exact results (if you evaluate the polynomial with them you will get a genuine zero), the one with RDF gives you approximate result.

EDIT : In your particular case, it seems that Sage is able to find the roots of the symbolic polynomial:

sage: solve(P(x), x)
[x == -1/2*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + 1/9*(8*I*sqrt(3) - 8)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 1/3, x == -1/2*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(-I*sqrt(3) + 1) + 1/9*(-8*I*sqrt(3) - 8)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 1/3, x == (1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 16/9/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 1/3, x == -1]


However, the result is far from being trustable and i would not recommend to use it in general, for example there will be situations where solve ...

more

@tmonteil The problem is that (even in your last segment) Sage seems to be finding complex roots! This is not possible! And is there a way to get this polynomial - the one I quoted in my question? (...the P that you have has constants as coeficients - what do you mean by "symbolic coefficients" ? ...)

( 2015-01-21 18:41:52 +0200 )edit

This is not true: when we dealt with genuine polynomials over the rings ZZ, QQ, QQbar, RDF, we only got real roots. The only place where we got apparently non-real complex roots is the last edit which was explicitly not recommended. Indeed, you got symbolic expressions, a framework where it is hard to know whether an element is zero. In particular, Sage is not able to deduce that the first root is real:

sage: e = -1/2*(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*(I*sqrt(3) + 1) + 1/9*(8*I*sqrt(3) - 8)/(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) + 1/3
sage: e.imag_part()
-1/2*sqrt(3)*abs(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)*cos(1/3*arctan2(1/9*sqrt(101)*sqrt(3), 37/27)) - 1/2*abs(1/9*I*sqrt ...
(more)
( 2015-01-21 23:25:57 +0200 )edit