#20019: Wrong approximate real roots of cubic equation and matrix eigenvalues
-----------------------------+-------------------------------------
Reporter: dlaurain | Owner:
Type: PLEASE CHANGE | Status: new
Priority: minor | Milestone: sage-7.1
Component: numerical | Keywords: eigenvalues cubic
Merged in: | Authors:
Reviewers: | Report Upstream: N/A
Work issues: | Branch:
Commit: | Dependencies:
Stopgaps: |
-----------------------------+-------------------------------------
I compute real roots of cubic equation in order to get intersection of
conics (matrix representation of conics C1,C2 and trying to find x such as
C1 + xC2 is a degenerate conic in the pencil of conics Cx, leading
computation eigenvalues of a 3x3 symmetric matrix).
The cubic equation is the characteristic polynomial of the matrix
representing the conic.
I noticed very annoying points : real elements of the matrix are
approximated to rational or square root of rational, then (via Maxima) a
usual direct formula is computed for roots involving terms with exponent
1/3, and results are not pure real '''but with a small imaginary part
(xxxe-15*I or xxxe-16*I)''' ...failing the is_real is_test.
Algorithm is wrong somewhere (approximate results must be reals), maybe in
the one third-root extracting.
I will bypass this myself for the special case cubic root solving (one
method given in wikipedia).
I vote for a removal of the imaginary part and an option for eigenvalues()
in order NOT to convert into the rational-squared rational numbers (which
is BIG computer-time consuming)
Here is the code and result (copy to sage cloud sagews cell and run it) :
{{{
m = matrix(SR,3,3, [ [ -0.5898846112288165, -1.2780028056044432,
0.19570118737659], [ -0.6573758889893906, 0.7138094643959463,
0.5958093145901924], [-0.15581144445167106,
-1.128948696332049,0.02254152681632878] ])
print "matrix m = ",m
e = m.eigenvalues()
print "e = ",e
for i in xrange(3):
print "e[",i,"] = ",e[i].n()
k = e[i].n()
if k.is_real():
print "real"
print "ok"
}}}
and the result:
{{{
matrix m = [ -0.589884611228817 -1.278002805604443 0.195701187376590]
[ -0.657375888989391 0.713809464395946 0.595809314590192]
[-0.1558114444516711 -1.12894869633205 0.0225415268163288]
e =
[-1/2*(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)*(I*sqrt(3)
+ 1) -
1/31319563713233492557164909561465418607943309241358730913107628212618325017905121528*(-2935792338917283029273723965573487435778349039987383711409958811428486614200489691*I*sqrt(3)
+
2935792338917283029273723965573487435778349039987383711409958811428486614200489691)/(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)
+ 1522704526641019831930/31188820126768548875109,
-1/2*(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)*(-I*sqrt(3)
+ 1) -
1/31319563713233492557164909561465418607943309241358730913107628212618325017905121528*(2935792338917283029273723965573487435778349039987383711409958811428486614200489691*I*sqrt(3)
+
2935792338917283029273723965573487435778349039987383711409958811428486614200489691)/(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)
+ 1522704526641019831930/31188820126768548875109,
(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)
+
2935792338917283029273723965573487435778349039987383711409958811428486614200489691/15659781856616746278582454780732709303971654620679365456553814106309162508952560764/(1/1343995347401075311202748407859241110077609781036218957900022883378876944992630605072050285851332957238651732936623207509906080896073076981060324480155982437131228948824*I*sqrt(1524851746879596218535094385454836889563018461734179933863331424825683633061404780241575261901618255004533012191704903371946279282869850682188819701138931149893341570753485037882769059356757126038265431530617860090414719193317893)*sqrt(4075786085992596757393163640405535653407766182130452275086346789770570836250529766467242936781545895842813)
-
15903655260910485632495073834892766895333845594295595601980518211356549746527605047218514971863402134557858900217697755/283437780734931771868343391472037014998470365820422957039812196088793940426823102967471008021664106311557338766483164812)^(1/3)
+ 1522704526641019831930/31188820126768548875109]
e[ 0 ] = 0.266756311146258 - 1.88737914186277e-15*I
e[ 1 ] = -0.785953947604597 + 5.55111512312578e-16*I
e[ 2 ] = 0.665664016441800 + 1.33226762955019e-15*I
ok
}}}
--
Ticket URL: <http://trac.sagemath.org/ticket/20019>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/d/optout.
