Cubic equations

[Rainer Dorsch]

Hello,

when solving a cubic equation

fk:  95.9271531456121*x - 1.7e-18/(x - 6.0e-7)**2

I get a strange result:

[6.10000000000001e-8 - 0.e-30*I, 3.85627081928849e-7 + 0.e-27*I,
7.53372918071151e-7 - 0.e-28*I]

What does the 0.e-30*I mean? I do not think that there should be an imaginary
part...

Even worse, if I decleare x as real, I do get no solution at all.


Here is some reference code

import sympy as sp
sp.init_printing()


k,x=sp.symbols('k,x')
#k,x=sp.symbols('k,x',real=True)

f=k*x - 1.7e-18/(x - 6e-7)**2

x_test=6.1e-8
print("x_test set: ",x_test)
fx=f.subs(x,x_test)
print("fx: ",fx)
k0=sp.solvers.solve(fx,k)[0]
print("k: ",k0)

print("x_test:")
fk=sp.simplify(f.subs(k,k0))
print("fk: ",fk)

p1=sp.plotting.plot(fk,(x,0,1e-7))

x_test=sp.solvers.solve(fk,x)
print(x_test)

and the output

*** example end ***

(sympy) rd@h370-wlan:~$ python3 online2.py  
x_test set:  6.1e-08
fx:  6.1e-8*k - 5.85155634188234e-6
/home/rd/virtualenv/sympy/lib/python3.7/site-packages/sympy/__init__.py:676:
SymPyDeprecationWarning:  

importing sympy.solvers.solvers with 'from sympy import *' has been
deprecated since SymPy 1.6. Use import sympy.solvers.solvers instead.
See https://github.com/sympy/sympy/issues/18245 for more info.

 deprecated_since_version="1.6").warn()
k:  95.9271531456121
x_test:
fk:  95.9271531456121*x - 1.7e-18/(x - 6.0e-7)**2
2.8e-06 |                                                     ..
| ...  
| ...      
| ...        
| ..            
| ...              
| ...                
| ...                    
| ..                      
| ...                        
-9.5e-0 |-------------------------...---------------------------
| ..                              
| ...                                
| ...                                    
| ..                                      
| ...                                        
| ..                                            
| ...                                              
| ...                                                
| ..                                                    
-4.7e-0 |_______________________________________________________
        0                          5e-08                      1e-07
[6.10000000000001e-8 - 0.e-30*I, 3.85627081928849e-7 + 0.e-27*I,
7.53372918071151e-7 - 0.e-28*I]
(sympy) rd@h370-wlan:~$

***** Output end *****


Any hint is welcome
Thanks
Rainer

[Aaron Meurer]

The very small imaginary part is due to floating point issues. You can
ignore it, or remove it with solution.evalf(chop=True).

The reason this shows up is because SymPy uses the cubic equation,
which requires imaginary numbers to express, even when the full root
is real (https://en.wikipedia.org/wiki/Casus_irreducibilis). This
means that the value is mathematically real, but due to floating point
inaccuracies, the floating point value ends up with a small imaginary
part.

Setting real=True causes SymPy to filter out non-real solutions. It
isn't smart enough to know that 6.10000000000001e-8 - 0.e-30*I, should
be real because it is only doing a simple filtering on the roots
without any consideration of the original equation.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHXQ_BJj-dyPUiNZRi0bVSFPK1Cmi-VVnC5sDYc35GjwZg-q2g%40mail.gmail.com.

[David Bailey]

On 28/05/2020 22:13, Rainer Dorsch wrote:

Hello,

when solving a cubic equation

fk:  95.9271531456121*x - 1.7e-18/(x - 6.0e-7)**2

I get a strange result:

[6.10000000000001e-8 - 0.e-30*I, 3.85627081928849e-7 + 0.e-27*I,
7.53372918071151e-7 - 0.e-28*I]

What does the 0.e-30*I mean? I do not think that there should be an imaginary
part...

Even worse, if I decleare x as real, I do get no solution at all.

In addition to Aaron's comments I would like to point  out that if you have a non-linear equation which is either intrinsically inexact because one or more of the coefficients is a floating point number (as in this case) or where the intention is simply to convert the answer as a floating point number, it is far better to use nsolve in the real domain:

nsolve( 95.9271531456121*x - 1.7e-18/(x - 6.0e-7)**2,x,(-50,50))

The result is generally more accurate, spurious imaginary solutions don't appear, and the calculation is probably more efficient (only relevant if the intention is to do a lot of such calculations).

David

[Oscar Benjamin]

There are also specific numerical algorithms for finding approximate
roots of univariate polynomials the can be used through the roots
function or the nroots method of Poly:

>>> from sympy import *
>>> x = Symbol('x')
>>> eq = 95.9271531456121*x - 1.7e-18/(x - 6.0e-7)**2
>>> eq2 = eq.as_numer_denom()[0]
>>> eq2
95.9271531456121*x*(x - 6.0e-7)**2 - 1.7e-18
>>> roots(eq2)
{6.10000000000001e-8: 1, 3.85627081928849e-7: 1, 7.53372918071150e-7: 1}
>>> p = Poly(eq2, x)
>>> p
Poly(95.9271531456121*x**3 - 0.000115112583774735*x**2 +
3.45337751324203e-11*x - 1.7e-18, x, domain='RR')
>>> p.nroots()
[6.10000000000001e-8, 3.85627081928849e-7, 7.53372918071150e-7]

Because these are tailored to polynomials they can have explicit
knowledge about which roots are real and which are non-real:

>>> Poly(x**3 + x**2 + x + 1, x).nroots()
[-1.00000000000000, -1.0*I, 1.0*I]

--
Oscar