CODEWITHSUNDEEP
numpy
user2820579
user2820579

Reputation: 3451

Strange roots `using numpy.roots`

Is there something wrong in the evaluation of the polinomial (1-alpha*z)**9 using numpy? For alpha=3/sqrt(2) my list of coefficients is given in the array

psi_t0 =  [1.0, -19.0919, 162.0, -801.859, 2551.5, -5412.55, 7654.5, -6958.99, 3690.56, -869.874]

According to numpy documentation, I have to invert this array in order to compute the zeros, i.e.

psi_t0 = psi_t0[::-1]

Thus giving

a = np.roots(psi_t0)
[0.62765842+0.06979364j 0.62765842-0.06979364j 0.52672941+0.14448097j  0.52672941-0.14448097j 0.42775926+0.13031547j 0.42775926-0.13031547j  0.36690056+0.07504044j 0.36690056-0.07504044j 0.34454214+0.j]

which is completely crap since the roots must be all equal to sqrt(2)/3.

Upvotes: 1

Views: 190

Answers (1)

Paul Panzer
Paul Panzer

Reputation: 53029

As you take the 9th power you'll find that you create a very "wide" zero, indeed, if you step eps away from the true zero and evaluate you'll get something of O(eps^9). In view of that numerical inaccuracies are all but expected.

>>> np.set_printoptions(4)
>>> print(C)
[-8.6987e+02  3.6906e+03 -6.9590e+03  7.6545e+03 -5.4125e+03  2.5515e+03
 -8.0186e+02  1.6200e+02 -1.9092e+01  1.0000e+00]
>>> np.roots(C)
array([0.4881+0.0062j, 0.4881-0.0062j, 0.4801+0.0154j, 0.4801-0.0154j,
       0.4681+0.0172j, 0.4681-0.0172j, 0.458 +0.011j , 0.458 -0.011j ,
       0.4541+0.j    ])
>>> np.polyval(C,_)
array([1.4622e-13+6.6475e-15j, 1.4622e-13-6.6475e-15j,
       1.2612e-13+1.5363e-14j, 1.2612e-13-1.5363e-14j,
       1.0270e-13+1.3600e-14j, 1.0270e-13-1.3600e-14j,
       1.1346e-13+9.7179e-15j, 1.1346e-13-9.7179e-15j,
       1.0936e-13+0.0000e+00j])

As you can see the roots numpy returns are "good" in that the polynomial evaluates to something pretty close to zero at these points.

Upvotes: 2

Related Questions