ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 21 Jan 2015 16:25:57 -0600How to stop Sage from finding erroneous complex roots?http://ask.sagemath.org/question/25587/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?
Tue, 20 Jan 2015 23:06:19 -0600http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/Answer by tmonteil for <p>Consider the matrix, </p>
<pre><code>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])
</code></pre>
<p>where w is the m^th of the k^th roots of unity, <code>w = exp((2*pi*I*m )/k</code>
for some k a positive integer, and 1 <= m <= (k-1)
and 1<= a,b,c <= (k-1) </p>
<p>Then the characteristic polynomial of the above matrix is,
<code>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) )</code></p>
<ul>
<li>Is there a way to get sage to be able to calculate the above characteristic polynomial? </li>
</ul>
<p>Now I try getting roots of the above by doing,</p>
<pre><code>g(x)=real_part(p(x)).simplify()
g.solve(x)
</code></pre>
<ul>
<li><p>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! </p>
<p>How to get across this trouble? </p></li>
</ul>
http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?answer=25589#post-id-25589The 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` will silently gives you only a (possibly empty) subset of the solutions
Wed, 21 Jan 2015 03:03:35 -0600http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?answer=25589#post-id-25589Comment by tmonteil for <div class="snippet"><p>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:</p>
<pre><code>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
</code></pre>
<p>You can compute the characteristic polynomial:</p>
<pre><code>sage: P = A.characteristic_polynomial()
sage: P
x^4 - 6*x^2 - 6*x - 1
</code></pre>
<p>But it has symbolic coefficients, so it is not able to find its roots:</p>
<pre><code>sage: P.parent()
Univariate Polynomial Ring in x over Symbolic Ring
sage: P.roots()
NotImplementedError:
</code></pre>
<p>What you should do is to specify a better ring for the coefficients. You will see that the roots depend on the ring:</p>
<p>There is one integer root with multiplicity 1:</p>
<pre><code>sage: Q = P.change_ring(ZZ)
sage: Q.roots()
[(-1, 1)]
</code></pre>
<p>There is one rational root with multiplicity 1:</p>
<pre><code>sage: Q = P.change_ring(QQ)
sage: Q.roots()
[(-1, 1)]
</code></pre>
<p>They are 4 algebraic roots, each with multipicity 1.</p>
<pre><code>sage: Q = P.change_ring(QQbar)
sage: Q.roots()
[(-1.655442381549831?, 1),
(-1, 1),
(-0.2107558809591918?, 1),
(2.866198262509023?, 1)]
</code></pre>
<p>You can also work in inexact rings, to get floating point approximations, for example in the <code>Real Double Field</code>:</p>
<pre><code>sage: Q = P.change_ring(RDF)
sage: Q.roots()
[(-1.6554423815498323, 1),
(-0.9999999999999978, 1),
(-0.21075588095919173, 1),
(2.866198262509024, 1)]
</code></pre>
<p>Note that there is a huge difference between the last two blocks, the one with <code>QQbar</code> gives you exact results (if you evaluate the polynomial with them you will get a genuine zero), the one with <code>RDF</code> gives you approximate result.</p>
<p><strong>EDIT</strong> : In your particular case, it seems that Sage is able to find the roots of the symbolic polynomial:</p>
<pre><code>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]
</code></pre>
<p>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 <code>solve ...</code><span class="expander"> <a>(more)</a></span></p></div>http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?comment=25599#post-id-25599This 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(101)*sqrt(3) + 37/27)^(1/3)*sin(1/3*arctan2(1/9*sqrt(101)*sqrt(3), 37/27)) + 8/9*sqrt(3)*cos(-1/3*arctan2(1/9*sqrt(101)*sqrt(3), 37/27))/abs(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3) - 8/9*sin(-1/3*arctan2(1/9*sqrt(101)*sqrt(3), 37/27))/abs(1/9*I*sqrt(101)*sqrt(3) + 37/27)^(1/3)
sage: bool(e.imag_part() == 0)
True
The expression of the imaginary part of `e` is so huge that Sage is not able to decide whether it is zero. However, you can find a numerical approximation that tends to prove that this expression has zero imaginary part:
sage: RDF(e.imag_part())
-2.498001805406602e-16Wed, 21 Jan 2015 16:25:57 -0600http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?comment=25599#post-id-25599Comment by Phoenix for <div class="snippet"><p>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:</p>
<pre><code>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
</code></pre>
<p>You can compute the characteristic polynomial:</p>
<pre><code>sage: P = A.characteristic_polynomial()
sage: P
x^4 - 6*x^2 - 6*x - 1
</code></pre>
<p>But it has symbolic coefficients, so it is not able to find its roots:</p>
<pre><code>sage: P.parent()
Univariate Polynomial Ring in x over Symbolic Ring
sage: P.roots()
NotImplementedError:
</code></pre>
<p>What you should do is to specify a better ring for the coefficients. You will see that the roots depend on the ring:</p>
<p>There is one integer root with multiplicity 1:</p>
<pre><code>sage: Q = P.change_ring(ZZ)
sage: Q.roots()
[(-1, 1)]
</code></pre>
<p>There is one rational root with multiplicity 1:</p>
<pre><code>sage: Q = P.change_ring(QQ)
sage: Q.roots()
[(-1, 1)]
</code></pre>
<p>They are 4 algebraic roots, each with multipicity 1.</p>
<pre><code>sage: Q = P.change_ring(QQbar)
sage: Q.roots()
[(-1.655442381549831?, 1),
(-1, 1),
(-0.2107558809591918?, 1),
(2.866198262509023?, 1)]
</code></pre>
<p>You can also work in inexact rings, to get floating point approximations, for example in the <code>Real Double Field</code>:</p>
<pre><code>sage: Q = P.change_ring(RDF)
sage: Q.roots()
[(-1.6554423815498323, 1),
(-0.9999999999999978, 1),
(-0.21075588095919173, 1),
(2.866198262509024, 1)]
</code></pre>
<p>Note that there is a huge difference between the last two blocks, the one with <code>QQbar</code> gives you exact results (if you evaluate the polynomial with them you will get a genuine zero), the one with <code>RDF</code> gives you approximate result.</p>
<p><strong>EDIT</strong> : In your particular case, it seems that Sage is able to find the roots of the symbolic polynomial:</p>
<pre><code>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]
</code></pre>
<p>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 <code>solve ...</code><span class="expander"> <a>(more)</a></span></p></div>http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?comment=25595#post-id-25595@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" ? ...)Wed, 21 Jan 2015 11:41:52 -0600http://ask.sagemath.org/question/25587/how-to-stop-sage-from-finding-erroneous-complex-roots/?comment=25595#post-id-25595