[sage-support] Re: Spurious numerical solutions of polynomial equation
Curious. The polynomial does not seem particularly ill-conditioned, in the
sense that its discriminat is roughly what one might expect (unlike, say,
Wilkinson's). Maple gives 50 roots with |f(x)|10^{-12}, and one with f(x)=,
i.e. much larger.
On Thursday, 12 December 2013 15:35:53 UTC, AWWQUB wrote:
Consider the equation
*(I*x^51+sum(x^k,k,0,50))==0*
Try to solve it numerically using
*solve([(I*x^51+sum(x^k,k,0,50))==0,x==x],x,solution_dict=True)*
and you obtain 51 solutions of which 50 have modulus approximately 1 and
the other is close to 1+I. Substituting back gives residuals of around at
most 10^-6. That looks fine and Mathematica gives similar solutions. Now
substitute x-1-I for x, so that one of the solutions should now be close to
zero. Sage now gives solutions with real and imaginary parts both between
-1 and +4, but none are particularly close to zero.The residuals now range
up to 10^20. Mathematica gives completely different answers but which are
also wrong, namely all 51 are approximately 1.
Any ideas?
Tony Wickstead
--
You received this message because you are subscribed to the Google Groups
sage-support group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/groups/opt_out.
[sage-support] Re: Spurious numerical solutions of polynomial equation
Apparently you need about 65 bits of accuracy, presumably maxima uses
doubles (=53 bits). Compare:
sage: eq = (I*x^51+sum(x^k,k,0,50)).subs(x=x-1-I)
sage: var('x, k')
(x, k)
sage: eq = (I*x^51+sum(x^k,k,0,50)).subs(x=x-1-I)
sage: sol = eq.polynomial(ComplexField(65)).roots()
sage: [abs(eq.subs(x=root[0])) for root in sol]
[6.549903740982083181e-19,
3.535528103523514460e-19,
1.674320837445456618e-19,
2.471478998132169129e-19,
6.778825368452881617e-19,
1.898563463306179569e-19,
2.315856855599388456e-19,
6.564587311639428261e-19,
2.032879073410320814e-19,
8.506842597138909958e-19,
2.751697050590121811e-19,
3.836227086515876689e-19,
9.128610475519328751e-19,
6.517906145587456560e-19,
3.388131789017201356e-19,
4.824490284086275634e-19,
2.710505431213761085e-19,
5.826404650794643473e-19,
5.762613592002297537e-19,
3.095091306081366359e-19,
8.243680437578442218e-19,
4.85449649890104e-19,
2.498962532480909844e-19,
5.354380611373575798e-19,
6.465564001006996412e-19,
3.615352538357994336e-19,
8.292552051274916129e-19,
1.066116732848065361e-18,
9.482896068074615061e-19,
9.157989237069755039e-19,
3.794707603699265519e-19,
9.953556208659440424e-19,
1.148691607350736122e-18,
1.655905868183726679e-18,
1.540940812835139611e-18,
3.632102562304021384e-19,
1.379498152005247844e-18,
9.199760500100244489e-19,
2.275896670392357755e-18,
2.091835968936370797e-18,
1.894447462943825357e-18,
1.012728404406850756e-18,
2.548451743371480535e-18,
3.297472175031376887e-18,
4.882597681423819153e-18,
5.704181241051244702e-18,
6.093334370826141925e-18,
4.233944833136100806e-18,
9.844770653768304268e-18,
3.233122873428597620e-18,
1.601494307857450561e-14]
On Thursday, December 12, 2013 3:35:53 PM UTC, AWWQUB wrote:
Consider the equation
*(I*x^51+sum(x^k,k,0,50))==0*
Try to solve it numerically using
*solve([(I*x^51+sum(x^k,k,0,50))==0,x==x],x,solution_dict=True)*
and you obtain 51 solutions of which 50 have modulus approximately 1 and
the other is close to 1+I. Substituting back gives residuals of around at
most 10^-6. That looks fine and Mathematica gives similar solutions. Now
substitute x-1-I for x, so that one of the solutions should now be close to
zero. Sage now gives solutions with real and imaginary parts both between
-1 and +4, but none are particularly close to zero.The residuals now range
up to 10^20. Mathematica gives completely different answers but which are
also wrong, namely all 51 are approximately 1.
Any ideas?
Tony Wickstead
--
You received this message because you are subscribed to the Google Groups
sage-support group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/groups/opt_out.
[sage-support] Re: Spurious numerical solutions of polynomial equation
Sorry - I was looking at the original *(I*x^51+sum(x^k,k,0,50))==0.* Note
that if you make the x-x-1 substitution, the polynomial now has
coefficient as large as 10^18. The discriminant is unchanged, but the value
you would expect it to be, given the size of the coefficient, is
roughly10^800 times greater, so the polynomial is now amazingly
ill-conditioned. Note that if you evaluate the translated polynomial in
machine floating-point, and then translate back, you get coefficients of
size 10^15. This translation is not feasible in floating-point!
On Saturday, 14 December 2013 17:37:16 UTC, JamesHDavenport wrote:
Curious. The polynomial does not seem particularly ill-conditioned, in the
sense that its discriminant is roughly what one might expect (unlike, say,
Wilkinson's). Maple gives 50 roots with |f(x)|10^{-12}, and one with f(x)=,
i.e. much larger.
On Thursday, 12 December 2013 15:35:53 UTC, AWWQUB wrote:
Consider the equation
*(I*x^51+sum(x^k,k,0,50))==0*
Try to solve it numerically using
*solve([(I*x^51+sum(x^k,k,0,50))==0,x==x],x,solution_dict=True)*
and you obtain 51 solutions of which 50 have modulus approximately 1 and
the other is close to 1+I. Substituting back gives residuals of around at
most 10^-6. That looks fine and Mathematica gives similar solutions. Now
substitute x-1-I for x, so that one of the solutions should now be close to
zero. Sage now gives solutions with real and imaginary parts both between
-1 and +4, but none are particularly close to zero.The residuals now range
up to 10^20. Mathematica gives completely different answers but which are
also wrong, namely all 51 are approximately 1.
Any ideas?
Tony Wickstead
--
You received this message because you are subscribed to the Google Groups
sage-support group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/groups/opt_out.
