Re: Re: Beware of NSolve - nastier example

*To*: mathgroup at smc.vnet.net*Subject*: [mg50372] Re: [mg50346] Re: Beware of NSolve - nastier example*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Tue, 31 Aug 2004 06:28:50 -0400 (EDT)*References*: <200408200858.EAA12533@smc.vnet.net> <cg6srb$odf$1@smc.vnet.net> <200408280837.EAA19074@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Carlos Felippa wrote: > Daniel Lichtblau <danl at wolfram.com> wrote in message news:<cg6srb$odf$1 at smc.vnet.net>... > < ... > > >>On the flip side, the heuristics that decide when a root actually >>satisfies the equation can be a bit lax. Hence parasites can wind up in >>the returned solution set. One way to influence this is to rescale the >>equation e.g. multiplying by some large number. >> >>I'll look into tightening some of the decisions as to what constitutes a >>sufficiently small residual. The near-double-root issue goes with the >>territory and is in no way incorrect behavior. >> >>What is NR, by the way? >> >> >>Daniel Lichtblau >>Wolfram Research > > > Newton and Raphson from Olde England. Joseph Raphson was > rumored to be Newton's programmer. > > BTW your suggestion to scale the equation does yield some > improvements -- at least no wrong roots appear: > > f=10^8*(5/432-11/(27*Sqrt[70]*Sqrt[19-1890*x])+x/(2*Sqrt[38/35-(108+1/10000000)*x])); > Print[N[Solve[f==0,x]]]; (* gives 3 roots *) > Print[NSolve[f,x,16]]; (* 3 correct roots *) > Print[NSolve[f,x,21]]; (* 2 roots, missed 1 *) > Print[NSolve[f,x,24]]; (* 2 roots, missed 1 *) > Print[NSolve[f,x,28]]; (* 2 roots, missed 1 *) > Print[NSolve[f,x,32]]; (* 2 roots, missed 1 *) > Print[NSolve[f,x,64]]; (* 2 roots, missed 1 *) > Print[NSolve[f,x,128]]; (* 3 correct roots *) As for that ill conditioned root, let me say a bit more. There are four important distinct values that agree to single precision. They are the actual root closest to 1/100, a parasite root, and the roots of the denominator quantities Sqrt[19-1890*x] and Sqrt[38/35-(108+1/10000000)*x]. The parasite should not cause trouble but the denominator singularities will make it effectively impossible to get a solution at single precision except perhaps by luck. Here are some experiments I tried along the lines of iterative solvers. f = 10^8*(5/432-11/(27*Sqrt[70]*Sqrt[19-1890*x])+x/ (2*Sqrt[38/35-(108+1/10000000)*x])); First do FindRoot starting very close to the actual root. In[52]:= InputForm[FindRoot[f==0, {x, 0.01005291004146810952497038}]] FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. Out[52]//InputForm= {x -> 0.01005291004146811} It stayed there, which is good. Now try nearer to the parasite root. In[53]:= InputForm[FindRoot[f==0, {x, 0.01005291004146847618464402}]] FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. Out[53]//InputForm= {x -> 0.01005291004146811} It moved nearer to the actual root, which is good. Not a surprise though because the actual and parasite roots are far closer together than they are to the denominator singularities. This time we start nearer to one of those singularities. In[55]:= InputForm[FindRoot[f==0, {x, 0.01005291004360180287}]] 1 Power::infy: Infinite expression -------- encountered. Sqrt[0.] FindRoot::nlnum: The function value {ComplexInfinity} is not a list of numbers with dimensions {1} at {x} = {0.0100529}. Out[55]//InputForm= FindRoot[f == 0, {x, 0.01005291004360180287`18.002291796377463}] We are still near the singularity and obviously the blowup caused trouble (as well it ought). Now let's look at a few attempts with a single precision starting point. In[56]:= InputForm[FindRoot[f==0, {x, 0.01}, Method->Newton]] Out[56]//InputForm= {x -> 0.0028165928963646216} In[57]:= InputForm[FindRoot[f==0, {x, 0.01,0.0101}, Method->Secant]] Out[57]//InputForm= {x -> 0.0028165928963646133 + 1.8983848342175145*^-19*I} In[58]:= InputForm[FindMinimum[Abs[f], {x, 0.01}, Method->Newton]] FindMinimum::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. Out[58]//InputForm= {0.18022000758052825, {x -> 0.0028165851990656406}} In[59]:= InputForm[FindMinimum[Abs[f], {x, 0.01}, Method->QuasiNewton]] FindMinimum::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. Out[59]//InputForm= {4.83554168928535*^-9, {x -> 0.0028165928963648224}} Reliably we seem to get to a different root (albeit with sometimes high residual, which indicates why scaling might be useful). Here I do simple Newton iterations keeping precision held constant. In[60]:= Block[ {$MinPrecision=20,$MaxPrecision=20}, x1 = N[1/100,20]; Do[Print[x1]; x1 = x1 - (f/D[f,x] /. x->x1), {10}] ] 0.010000000000000000000 0.0099073848257359387420 0.0096744232229347655343 0.0091485713560636477336 0.0081362079081019597712 0.0065730142014561307872 0.0047762592170251683323 0.0034097065684361748104 0.0028780802088755471897 0.0028172949584050694447 Same story. It's clear that this will not converge to the root closest to 1/100. Hence my claim that one root can be expected to cause extreme duress for iterative solvers at low precision. Daniel Lichtblau Wolfram Research

**References**:**Beware of NSolve - nastier example***From:*carlos@colorado.edu (Carlos Felippa)

**Re: Beware of NSolve - nastier example***From:*carlos@colorado.edu (Carlos Felippa)

**Re: Question: weight matrices and MonomialOrder for GroebnerBasis**

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

**3D Screensavers for Mathematica users**