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

*To*: mathgroup at smc.vnet.net*Subject*: [mg50390] Re: [mg50372] Re: [mg50346] Re: Beware of NSolve - nastier example*From*: DrBob <drbob at bigfoot.com>*Date*: Wed, 1 Sep 2004 01:49:35 -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>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

Oops, in my last paragraph I claimed convergence for a wider range of starting values than I actually get. (I had a typo when doing some of the tests.) 0.003`n to 0.04`n looks fine. At least it converges for SOME input, which is more than the raw Newton iteration can do. Of course, I used Solve to get the h function--the same Solve routine that solved the original problem already! h is expressible in radicals, though, where the original Solve result was not. Bobby On Tue, 31 Aug 2004 11:34:15 -0500, DrBob <drbob at bigfoot.com> wrote: > Here's a single-precision solution that's exactly right to machine precision. But look at the function value! > > f[x_] = 10^8*(5/ > 432 - 11/( > 27*Sqrt[70]*Sqrt[ > 19 - 1890*x]) + x/(2*Sqrt[38/35 - (108 + 1/10000000)*x])); > single = N[x /. Last@Solve[f@x == 0, x]] > f@single > > 0.01005291004146811 > 6767.298384602327 > > Is that the solver's fault? No, the Solver is dead on. This is just a single precision rendition of the exact answer. Here's a 50-digit answer: > > fifty = N[x /. Last@Solve[f@x == 0, x], 50] > f@fifty > > 0.010052910041468109524970376989018540981826996415272 > 0``29.97022428635945 > > Same solver, same answer, different precision. > > What about Newton iteration? Newton iterations work when the following function is a contraction: > > g[x_] := x - f'[x] > g'[fifty] > > -5.2651323530359985064929525726862151829251964942846117`39.90793170384201*^33 > > g needs a derivative strictly between -1 and 1, and that's QUITE a bit bigger. > > So there's no chance Newton iterations could work in any language, with any arithmetic model, any way you fix it. > > HOWEVER, if we replace the coefficient (x) of 1080000001 with y, and solve for y in terms of x, we do get a contraction near the root we're trying for. (I tried replacing the other mentions of x also, but they don't yield a contraction.) > > h[x_] = y /. First[Solve[100000000*(5/432 - 11/(27*Sqrt[70]* Sqrt[19 - 1890*x]) + > x/(2*Sqrt[38/35 - (1080000001*y)/10000000])) == 0, y]] > (4000000*(-610147 + 16720*Sqrt[70]*Sqrt[19 - 1890*x] + 31421250*x + 542959200* > x^2 - 54010152000*x^3))/(7560000007*(-32113 + 880*Sqrt[70]*Sqrt[19 - 1890*x] + 1653750*x)) > > h'[fifty] > 0.186486495352634630910534650649828034524831032248`35.538455152501946 > > A Plot suggests the derivative is consistently small, but single-precision inputs cause Plot to generate errors. There are one or more denominators in h[x] that are problematic. Iteration works, however. For instance: > > n = 10; > Block[{$MinPrecision = n, $MaxPrecision = n}, > Chop@FixedPoint[h, .005`5] > ] > f[%] > 0.0100529100413456552939875214`10.000000000000004 > 0``-9.93478116473507 > > Machine precision doesn't work at all, but fixed-precision works even if it's LESS than machine precision. I got convergence for a fairly wide range of initial starts, too, from 0.`5 to 1.0`5. (The initial precision could be anything at all, so long as it's specified.) > > Bobby > > On Tue, 31 Aug 2004 06:28:50 -0400 (EDT), Daniel Lichtblau <danl at wolfram.com> wrote: > >> 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 >> >> >> > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**newbie is looking for a customDistribution function**

**Inflight magazine puzzle**

**Re: newbie is looking for a customDistribution function**

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