2
$\begingroup$

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.

Thanks in advance 

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
Improve this question
$\endgroup$
2
  • $\begingroup$ I guess that machine precision algorithms are bound to fail with these high exponents and with this wide range of coefficients. $\endgroup$
    – Henrik Schumacher
    Feb 24, 2018 at 19:46
  • 3
    $\begingroup$ 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. $\endgroup$
    – Daniel Lichtblau
    Feb 24, 2018 at 19:57

1 Answer 1

Reset to default
2
$\begingroup$ 
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 *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Sandro
    Feb 25, 2018 at 12:39
  • $\begingroup$ @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. $\endgroup$
    – Bob Hanlon
    Feb 25, 2018 at 14:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.