Print real roots only in numpy

Question

I have something like this:

coefs = [28, -36, 50, -22]
print(numpy.roots(coefs))

Of course the result is:

[ 0.35770550+1.11792657j  0.35770550-1.11792657j  0.57030329+0.j ]

However, by using this method, how do I get it only to print the real roots if any (as floats)? Meaning just this for my example:

0.57030329

Answer 1

Do NOT use .iscomplex() or .isreal(), because roots() is a numerical algorithm, and it returns the numerical approximation of the actual roots of the polynomial. This can lead to spurious imaginary parts, that are interpreted by the above methods as solutions.

Example:

# create a polynomial with these real-valued roots:
p = numpy.poly([2,3,4,5,56,6,5,4,2,3,8,0,10])
# calculate the roots from the polynomial:
r = numpy.roots(p)
print(r) # real-valued roots, with spurious imaginary part
array([ 56.00000000 +0.00000000e+00j,  10.00000000 +0.00000000e+00j,
         8.00000000 +0.00000000e+00j,   6.00000000 +0.00000000e+00j,
         5.00009796 +0.00000000e+00j,   4.99990203 +0.00000000e+00j,
         4.00008066 +0.00000000e+00j,   3.99991935 +0.00000000e+00j,
         3.00000598 +0.00000000e+00j,   2.99999403 +0.00000000e+00j,
         2.00000000 +3.77612207e-06j,   2.00000000 -3.77612207e-06j,
         0.00000000 +0.00000000e+00j])
# using isreal() fails: many correct solutions are discarded
print(r[numpy.isreal(r)])
[ 56.00000000+0.j  10.00000000+0.j   8.00000000+0.j   6.00000000+0.j
   5.00009796+0.j   4.99990203+0.j   4.00008066+0.j   3.99991935+0.j
   3.00000598+0.j   2.99999403+0.j   0.00000000+0.j]

Use some threshold depending on your problem at hand instead. Moreover, since you're interested in the real roots, keep only the real part:

real_valued = r.real[abs(r.imag)<1e-5] # where I chose 1-e5 as a threshold
print(real_valued)