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
$\endgroup$
  • $\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 '18 at 19:46
  • 2
    $\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 '18 at 19:57
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 *)
$\endgroup$
  • $\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 '18 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 '18 at 14:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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