Re: Re: Beware of NSolve - nastier example
- To: mathgroup at smc.vnet.net
- Subject: [mg50426] Re: [mg50414] Re: Beware of NSolve - nastier example
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 3 Sep 2004 03:35:16 -0400 (EDT)
- References: <200408200858.EAA12533@smc.vnet.net> <cg6srb$odf$1@smc.vnet.net> <200408280837.EAA19074@smc.vnet.net> <200408311028.GAA18653@smc.vnet.net> <opsdluzdvziz9bcq@monster.cox-internet.com> <ch3olv$33$1@smc.vnet.net> <200409020834.EAA02117@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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
systems.
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
- References:
- Re: Beware of NSolve - nastier example
- From: carlos@colorado.edu (Carlos Felippa)
- Re: Beware of NSolve - nastier example