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: Re: newbie is looking**

**Re: Re: newbie is looking foracustomDistributionfunction**

**Re: Re: Beware of NSolve - nastier example**

**Inflight magazine puzzle**