Re: Possible bug in NSolve[equation, variable, precission]
- To: mathgroup at smc.vnet.net
- Subject: [mg74394] Re: [mg74299] Possible bug in NSolve[equation, variable, precission]
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Mon, 19 Mar 2007 22:06:23 -0500 (EST)
- References: <200703170710.CAA19095@smc.vnet.net>
Julian Aguirre wrote:
> Dear group,
>
> Mathematica 5.2 chokes solving numerically a polynomial equation.
>
> In[1] := $Version
> Out[1]= 5.2 for Mac OS X (64 bit) (June 20, 2005)
>
> In[2]:= poly=171142046150220198693105489-16023210221608713837587916
> x-2020825892011586434364754 x^2+190894692033395024364972
> x^3+6039743423966949379761 x^4-568929229651998950400
> x^5-470066550477520896 x^6+2821109907456 x^7;
>
> In[3]:= poly2=Expand[poly/9];
>
> In[4]:= NSolve[poly==0,x]
> Out[4]= {-1211.83, -13.0015, -13.0014, 11.923, 12.0809, 12.2509,
> 167826.}
>
> (* Up to this moment, everything is O.K. But *)
>
> In[5]:= NSolve[poly==0,x,20]
> Out[5]= $Aborted (* after a loooong time *)
>
> (* However, the following works as expected*)
>
> In[6]:= x/.NSolve[poly2==0,x,20]
> Out[6]= {-1211.8267955098487289, -13.001455891126, -13.001441554521,
> 11.92303189062617, 12.08089051352363, 12.25087466630727,
> 167826.26017849924816}
>
> Let me say that I have used Mathematica to solve thousands (probably
> millions) of equations like the one above. There must be some magic in
> the coefficients!
>
> Julian Aguirre
> University of the Basque Country
>
Yes, it has to do with the coefficients and some other features of the
polynomial. Adam Strzebonski and I investigated a similar set of
examples reported to us around two weeks ago. The user who brought it to
our attention also apparently does massive numbers of these in the
course of a large simulation, estimateds to be in the tens or hundreds
of thousands over a period of several hours. Over a period of many such
runs a few such examples were culled. So it does indeed seem to be a bit
of a rare phenomenon, though unfortunately not impossible to encounter.
What we know at present is this appears to involve a long-standing
problem in the numeric root-finding code used in the Mathematica kernel.
The trouble seems to correlate with near-multiplicity (but not with
actual multiplicity to within the requested bignum precision), along
with bad scaling. Your example also indicates it may have something to
do with the absolute magnitude of the coefficients, since rescaling by a
factor of 9 (the gcd of the coefficients) makes the problem disappear in
that case.
At this time we do not know in more detail what is behind the hang. One
thing to realize is that NSolve has dedicated special case handling for
a single univariate polynomial. Thus one workaround is to add an
equation in another variable e.g. y==x, to force NSolve to use a
different approach. Other possibilities include rescaling, or setting up
a companion matrix and finding the eigenvalues at whatever precision is
of interest (this uses Lapack methods rather than a Jenkins-Traub based
rootfinder).
Perhaps the main issue, for practical purposes, is how you can tell when
you will need a workaround. The best we have to offer is to use
TimeConstrained on your NSolve, and in cases where it aborts, go to plan b.
Daniel Lichtblau
Wolfram Research
- References:
- Possible bug in NSolve[equation, variable, precission]
- From: "Julian Aguirre" <julian.aguirre@ehu.es>
- Possible bug in NSolve[equation, variable, precission]