Sage returns complex solutions to equation with real solutions

[Oscar Gerardo Lazo Arjona]

I'm trying to find the solutions to solve this equation

sage: a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0
sage: b=solve(a==0,x)
sage: for i in b:
....: c=i.rhs()
....: print c.n()
....:
1.19783952189420 - 4.16333634234434e-17*I
0.998467807920659 + 1.38777878078145e-17*I
1.38386488335712 + 2.08166817117217e-17*I

But I get complex results instead of real ones. If I plot the function,
I can see that it definitely has 3 real solutions

plot(a,(x,0.9,1.5))

I realize that the imaginary parts of the results are minuscule, but
should they be there? In any case sage returns a symbolic expression,
and since I'm using floating point exponents I'd expect to get a
numerical aproximation to the solution.

Mathematica gets this right:

sage: mathematica_console()
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.

In[1]:= a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0;

In[2]:= Solve[a==0,x]

Out[2]= {{x -> 0.998468}, {x -> 1.19784}, {x -> 1.38386}}

thanks!

Oscar.

[Jeroen Demeyer]

On 2010-08-22 22:23, Oscar Gerardo Lazo Arjona wrote:
> I'm trying to find the solutions to solve this equation
>
> sage: a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0
> sage: b=solve(a==0,x)
> sage: for i in b:
> ....: c=i.rhs()
> ....: print c.n()
> ....:
> 1.19783952189420 - 4.16333634234434e-17*I
> 0.998467807920659 + 1.38777878078145e-17*I
> 1.38386488335712 + 2.08166817117217e-17*I
>
> But I get complex results instead of real ones. If I plot the function,
> I can see that it definitely has 3 real solutions

If you're only solving polynomial equations, this problem can be
overcome by defining a to be a polynomial:

sage: x = polygen(CC)
sage: a = 8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0
sage: parent(a)
Univariate Polynomial Ring in x over Complex Field with 53 bits of precision
sage: a.roots()
[(0.998467807920657, 1), (1.19783952189421, 1), (1.38386488335712, 1)]

Jeroen.

[Dima Pasechnik]

well, Mathematica probably does a counting of real roots before
solving
(and in any event an equation of odd degree with real coefficients is
guaranteed to get at least one
real solution).
It could also be that it just does a more aggressive rounding than
Sage.
In any event you can do this rounding, too, say

for i in b:
c=i.rhs()
if abs(c.n().imag())<1.e-15:
print c.n().real()
else: print c.n()

gives

1.19783952189420
0.998467807920659
1.38386488335712


On Aug 22, 10:23 pm, Oscar Gerardo Lazo Arjona

[cousteau]

On 22 ago, 22:23, Oscar Gerardo Lazo Arjona

<algebraicame...@gmail.com> wrote:
> I'm trying to find the solutions to solve this equation
>
> sage: a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0
> sage: b=solve(a==0,x)
> sage: for i in b:
> ....:     c=i.rhs()
> ....:     print c.n()
> ....:
> 1.19783952189420 - 4.16333634234434e-17*I
> 0.998467807920659 + 1.38777878078145e-17*I
> 1.38386488335712 + 2.08166817117217e-17*I

The problem with using real (floating point) numbers is that they're
probably solved numerically, using an algorithm that may be oriented
to finding also complex roots (the Laguerre method, for example). So
the first thing I'd suggest is entering the polynomial coefficients as
integer values:

sage: x = var('x'); A = 8594*x^3 - 30768*x^2 + 36399*x -14224
sage: B = solve(A, x)
sage: B[0]
x == -1/2*(I*sqrt(3) + 1)*(1/73856836*I*sqrt(2)*sqrt(5159598107) -
6552490/79340706073)^(1/3) - 1/73856836*(-457267*I*sqrt(3) + 457267)/
(1/73856836*I*sqrt(2)*sqrt(5159598107) - 6552490/79340706073)^(1/3) +
5128/4297

Ok, at least this might have an opportunity of getting simplified.
Let's try:

sage: B[0].full_simplify()
x == -1/17188*((I*sqrt(3) + 1)*(4297*I*sqrt(17)*sqrt(313)*sqrt(969667)
- 13104980*sqrt(2))^(2/3) -
10256*(4297*I*sqrt(17)*sqrt(313)*sqrt(969667) -
13104980*sqrt(2))^(1/3)*sqrt(2) - 457267*I*sqrt(3) + 457267)*sqrt(2)/
(4297*I*sqrt(17)*sqrt(313)*sqrt(969667) - 13104980*sqrt(2))^(1/3)

There doesn't seem to be a way to get rid of the I. Maybe Sage should
implement an algorithm for simplifying complex expressions and trying
to isolate real and imaginary parts, removing them if they're found to
be zero.


Back to floating point numbers, an alternative to solve() is
find_root(), which will attempt to find a numerical solution on a
given interval:

sage: a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0

sage: a.find_root(-5,5) # only finds a single root
1.383864883357113

[John Cremona]

It has been known since the 16th century that to solve a real cubic
equation with 3 real roots using the traditional formula involves the
use of complex numbers along the way! This was one reason for the
original acceptance of complex numbers.

In this case if you type a.roots? you will see that there are options,
and in fact a.roots(ring=RR) returns

sage: a.roots(ring=RR)

[(0.998467807920657, 1), (1.19783952189421, 1), (1.38386488335712, 1)]

or also

sage: a.roots(ring=RR, multiplicities=False)
[0.998467807920657, 1.19783952189421, 1.38386488335712]

However, I suggest that for many users who are not pure
mathematicians, having a different (or alternative?) name for the
parameter "ring" might be helpful.

John

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel+...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org
>

[Robert Dodier]

On Aug 22, 2:23 pm, Oscar Gerardo Lazo Arjona

<algebraicame...@gmail.com> wrote:

> sage: a=8594.0*x^3 - 30768.0 *x^2 + 36399.0 *x -14224.0
> sage: b=solve(a==0,x)
> sage: for i in b:
> ....:     c=i.rhs()
> ....:     print c.n()
> ....:
> 1.19783952189420 - 4.16333634234434e-17*I
> 0.998467807920659 + 1.38777878078145e-17*I
> 1.38386488335712 + 2.08166817117217e-17*I

IIRC Sage solve punts to Maxima solve, which tries to find
an exact solution. In Maxima I get

[x = (-sqrt(3)*%i/2-1/2)*(sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3)
+457267*(sqrt(3)*%i/2-1/2)/(36928418*(sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3))+5128/4297,
x = (sqrt(3)*%i/2-1/2)*(sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3)
+457267*(-sqrt(3)*%i/2-1/2)/(36928418*(sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3))+5128/4297,
x = (sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3)
+457267/(36928418*(sqrt(5159598107)*%i/
(18464209*2^(3/2))-6552490/79340706073)^(1/3))+5128/4297]

Now it's easy to see that these must all be real,
so essentially there must be an identity here (to make the imaginary
units go away) which is not being exploited by Maxima.

If you are specifically looking for numerical values,
probably there is a way to invoke a different algorithm in Sage.

FWIW

Robert Dodier

[rjf]

On Aug 23, 1:25 am, John Cremona <john.crem...@gmail.com> wrote:

> However, I suggest that for many users who are not pure
> mathematicians,  having a different (or alternative?) name for the
> parameter "ring" might be helpful.
>
> John
>

Maybe they shouldn't be using Sage if they don't know the term
"ring",
as well as the meaning of RR.

RJF

[Dr. David Kirkby]

If Sage wants to be a viable alternative to the 4 M's, then it needs to be able
to be usable by people who are not mathematicians. Despite having engineering
degrees I had never come across the term ring (in the mathematical sense),
before getting involved in Sage.

Dave

[kcrisman]

On Aug 23, 5:03 pm, "Dr. David Kirkby" <david.kir...@onetel.net>
wrote:

In the end, what is really going on here is that some very hard-
working individuals got a huge amount of basic functionality
implemented in a way where you don't need to know what a ring is, but
in ways that left lots of more user-intensive cases looking more
complicated than it wants to be. This just means much less hard-
working individuals (like me) eventually need to make this part of the
interface better.

But I think this *will* eventually happen. A lot of it will get
better when enough people understand Maxima well enough to make the
use of its flags etc. much better, which unfortunately has been on
hiatus for me for a while (as I remarked to Barton Willis off-list
earlier in one of these discussions). But it *will* happen.

- kcrisman

[rjf]

Maybe an alternative viewpoint (perhaps taken by Axiom, and maybe
Sage) is

(a) if you want to use a computer algebra system (especially one
accessed through Sage) effectively then
(b) you should know something about "modern algebra", at least so as
to be conversant with the notions
of ring, field, polynomial, ... and preliminary notions. This could
take about an hour of reading,
especially if you leave out proofs.

RJF

[Alex Ghitza]

And maybe we can even work some of this into the tutorial: some examples
of situations where you want to use the "ring" parameter together with a
note that "ring=number system" or even "ring=RR or CC or QQ".
Mathematicians reading this are not likely to be confused by it, and
non-mathematicians might not run away screaming.

People are likely to know what a polynomial *function* is, and we might
not need to make them aware that it's not the same thing as a
polynomial. But another bit that could go into the tutorial is some
good examples where you want to use polynomials instead of symbolics.

Best,
Alex

--
Alex Ghitza -- http://aghitza.org/
Lecturer in Mathematics -- The University of Melbourne -- Australia

[kcrisman]

Not at all. The viewpoint is rather related to a different thread:

(1) The people who develop most functionality originally are
mathematicians who do know about such things, and that was sufficient
for most purposes then, but

(2) Now we need to make those more user-accessible, and no one is
doing it for function X.

I agree this isn't ideal, but (as you know) open source software is
developed in a very different way, so we make the best of it rather
than giving up.

But there aren't really too many places where this comes in with
respect to "rings" etc.; it's really more of a system-wide thing where
there are options which you would know if you knew the underlying
software but aren't well-enough exposed at the Sage level (including
many Maxima ones Robert Marik has done an excellent job of uncovering
and making Pythonic syntax for).

- kcrisman