MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Beware of NSolve - nastier example

  • To: mathgroup at
  • Subject: [mg50426] Re: [mg50414] Re: Beware of NSolve - nastier example
  • From: Daniel Lichtblau <danl at>
  • Date: Fri, 3 Sep 2004 03:35:16 -0400 (EDT)
  • References: <> <cg6srb$odf$> <> <> <> <ch3olv$33$> <>
  • Sender: owner-wri-mathgroup at

Carlos Felippa wrote:
> [...]
> No problem whatsoever for conventional NR in 16-digits 
> (double precision) if you are sufficiently near a root.  Just run this
> NR iteration, which is at the level of the homeworks I assign to
> engineering sophomores:
> ClearAll[f,x];
> f=(5/432-11/(27*Sqrt[70]*Sqrt[19-1890*x])+x/(2*Sqrt[38/35-(108+1/10000000)*x]));
> fprime=D[f,x];
> Print[N[Solve[f==0,x]]//InputForm];  (* gives 3 correct roots *)
> SetPrecision[{xn,xnext,r,f,fprime,fp},16];
> xn=0.0100529100415;  (* the slippery root *)
> n=10;   (* actually n=3 is enough  to get limit accuracy *)
> For [i=1,i<=n,i++, r=N[f/.x->xn]; fp=N[fprime/.x->xn];
>      dx=-r/fp; xnext=xn+dx; 
>      Print[{i,r,fp,dx,xn,xnext}//InputForm]; xn=xnext];
> The derivative f' is large, O(10^13) whereas the residual goes down to
> O(10^(-4)). The interval where NR converges to that root is tiny, of
> O(10^(-10), but finite.  If one tries single precision (~ 6-7 digits)
> NR fails, as noted by Daniel Lichtblau,  since the convergence
> interval falls in the noise.
> What  is a bit surprising to me is that NSolve applied directly to f,
> (not as N[Solve][..]]) needs 128-digit working precision to resolve
> that root.

I can address this last part. It is a consequence of mathematical issues 
that arise in working with multiple or near multiple roots of algebraic 

First note that NSolve will augment with new variables (for radicals and 
denominators) and new equations (to record the algebraic dependencies 
between these variables). It arrives at a polynomial system, one that 
may have parasite roots it will later need to remove.

Second, I'll point out that for a root of multiplicity k, the eigencode 
that approximates the root values to n digits will typically have only 
n/k or so correct digits. So we must compromise in deciding whether a 
set of "nearby" values is actually a multiple root, because the methods 
used to find them will tend to smear out the multiple ones (to recover 
them, the usual rule is to average the values).

Recall that this particular algebraic system has a parasite root quite 
close to an actual root. NSolve needs to have enough precision to 
recognize that that root is distinct from the actual root. As they agree 
to 14 or so digits, and as there must be slop built into the code in 
order to recognize multiple roots at all, it requires far more than 14 
digits in order to make the determination that they are in fact distinct 
roots. Generally speaking, as the putative multiplicity is 2, roughly 30 
digits should suffice. But in practice we give it a larger multiplier. 
In this example 80 or so digits suffices.

So what goes wrong? First, so long as NSolve thinks it has a multiple 
root, the multiplicity reported would be incorrect. But, as I mentioned, 
the values will also be averaged as a correction step, and in the 
unfortunate situation where we do NOT have a multiple root that takes us 
away from the actual root. We then look at the residual. In this case it 
is not sufficiently small, at the precision utilized, for NSolve to 
retain it. So it is discarded as a parasite.

A possible way to address this in NSolve would be to first polish the 
root using e.g. FindRoot or FindMinimum. I have been reluctant to take 
that step because experimentation has shown it to be difficult to make 
reliable. One problem arises in your example. Suppose we take the 
approximated roots and use iterations to refine them. What if both 
approximations to the difficult root are in the attraction basin of the 
parasite instead of the actual root? What if both go to the actual root? 
Either way we get an incorrect result. It is only with sufficiently high 
precision approximations that we can reliably distinguish them, and it 
is not obvious to me how to get such precision from low precision 
approximations given the question of where they might converge in such 
refinement steps.

The overall picture is this. We have a global method for finding roots 
to algebraic systems. It works well but requires precision overestimates 
in order to handle parasite roots. In principle one might like it to 
apply local methods in the later stages. In practice I've been unable to 
make that work out.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Re: newbie is looking
  • Next by Date: Re: Re: newbie is looking foracustomDistributionfunction
  • Previous by thread: Re: Re: Beware of NSolve - nastier example
  • Next by thread: Inflight magazine puzzle