# problem in finding roots of a large polynomial

I have a problem with finding the roots of a huge polynomial with mathematica: it finds the roots but they are wrong!

I tried to use both "NSolve" as well as "Roots", but I always get the same wrong answer.

Can someone help to fix this? Below I attach the source.

P=1.05964*10^60 + 1.05858*10^57 s + 3.17279*10^59 s^2 + 3.06387*10^56 s^3 + 4.59167*10^58 s^4 + 4.28101*10^55 s^5 + 4.2773*10^57 s^6 + 3.84535*10^54 s^7 + 2.88161*10^56 s^8 + 2.49458*10^53 s^9 + 1.49555*10^55 s^10 + 1.24484*10^52 s^11 + 6.21943*10^53 s^12 + 4.96962*10^50 s^13 + 2.12826*10^52 s^14 + 1.62967*10^49 s^15 + 6.10696*10^50 s^16 + 4.47281*10^47 s^17 + 1.48993*10^49 s^18 + 1.04161*10^46 s^19 + 3.12282*10^47 s^20 + 2.07914*10^44 s^21 + 5.66691*10^45 s^22 + 3.58419*10^42 s^23 + 8.9553*10^43 s^24 + 5.36574*10^40 s^25 + 1.23757*10^42 s^26 + 7.00298*10^38 s^27 + 1.49987*10^40 s^28 + 7.98771*10^36 s^29 + 1.59677*10^38 s^30 + 7.97205*10^34 s^31 + 1.49409*10^36 s^32 + 6.96186*10^32 s^33 + 1.22805*10^34 s^34 + 5.31333*10^30 s^35 + 8.85215*10^31 s^36 + 3.53528*10^28 s^37 + 5.58004*10^29 s^38 + 2.04272*10^26 s^39 + 3.0631*10^27 s^40 + 1.01936*10^24 s^41 + 1.45581*10^25 s^42 + 4.36011*10^21 s^43 + 5.94408*10^22 s^44 + 1.58238*10^19 s^45 + 2.06352*10^20 s^46 + 4.8065*10^16 s^47 + 6.00697*10^17 s^48 + 1.19927*10^14 s^49 + 1.43889*10^15 s^50 + 2.39383*10^11 s^51 + 2.76176*10^12 s^52 + 3.67559*10^8 s^53 + 4.0836*10^9 s^54 + 407598. s^55 + 4.36684*10^6 s^56 + 290.57 s^57 + 3005.8 s^58 + 0.1 s^59 + 1. s^60

Sol = NSolve[P, s]

P /. Sol[[1, 1]]
2.45199*10^55 + 7.47092*10^54 I

• I guess that machine precision algorithms are bound to fail with these high exponents and with this wide range of coefficients. Feb 24, 2018 at 19:46
• Given the variance in scale of coefficients, it is not easy to see how one might accurately assess that a proposed root is incorrect. Testing with raised precision indicates that the roots found by NSolve actually are quite good, residuals notwithstanding. Feb 24, 2018 at 19:57

Clear[s]


First, use Rationalize to get exact numbers for the coefficients of the polynomial

P = 1.05964*10^60 + 1.05858*10^57 s + 3.17279*10^59 s^2 + 3.06387*10^56 s^3 +
4.59167*10^58 s^4 + 4.28101*10^55 s^5 + 4.2773*10^57 s^6 +
3.84535*10^54 s^7 + 2.88161*10^56 s^8 + 2.49458*10^53 s^9 +
1.49555*10^55 s^10 + 1.24484*10^52 s^11 + 6.21943*10^53 s^12 +
4.96962*10^50 s^13 + 2.12826*10^52 s^14 + 1.62967*10^49 s^15 +
6.10696*10^50 s^16 + 4.47281*10^47 s^17 + 1.48993*10^49 s^18 +
1.04161*10^46 s^19 + 3.12282*10^47 s^20 + 2.07914*10^44 s^21 +
5.66691*10^45 s^22 + 3.58419*10^42 s^23 + 8.9553*10^43 s^24 +
5.36574*10^40 s^25 + 1.23757*10^42 s^26 + 7.00298*10^38 s^27 +
1.49987*10^40 s^28 + 7.98771*10^36 s^29 + 1.59677*10^38 s^30 +
7.97205*10^34 s^31 + 1.49409*10^36 s^32 + 6.96186*10^32 s^33 +
1.22805*10^34 s^34 + 5.31333*10^30 s^35 + 8.85215*10^31 s^36 +
3.53528*10^28 s^37 + 5.58004*10^29 s^38 + 2.04272*10^26 s^39 +
3.0631*10^27 s^40 + 1.01936*10^24 s^41 + 1.45581*10^25 s^42 +
4.36011*10^21 s^43 + 5.94408*10^22 s^44 + 1.58238*10^19 s^45 +
2.06352*10^20 s^46 + 4.8065*10^16 s^47 + 6.00697*10^17 s^48 +
1.19927*10^14 s^49 + 1.43889*10^15 s^50 + 2.39383*10^11 s^51 +
2.76176*10^12 s^52 + 3.67559*10^8 s^53 + 4.0836*10^9 s^54 +
407598. s^55 + 4.36684*10^6 s^56 + 290.57 s^57 + 3005.8 s^58 + 0.1 s^59 +
1. s^60 // Rationalize[#, 0] &;


Secondly, use high precision in NSolve by using the option WorkingPrecision

sol = NSolve[P, s, WorkingPrecision -> 100];


Checking the results

P /. sol[[1 ;; 5]]

(* {0.*10^-23 + 0.*10^-23 I, 0.*10^-23 + 0.*10^-23 I, 0.*10^-21 + 0.*10^-21 I,
0.*10^-21 + 0.*10^-21 I, 0.*10^-26 + 0.*10^-26 I} *)

Max[Abs[P /. sol]]

(* 0.*10^-17 *)

• Many thanks! I did as you wrote and I am getting correct results. Just a curiosity about the command rationalize: what is the role of the "&" at the end of it? Feb 25, 2018 at 12:39
• @Sandro - It is a pure Function (&). Since the two argument form is needed, //Rationalize[#, 0]& is needed rather than the simpler //Rationalize used for a single argument. Feb 25, 2018 at 14:11