Re: Re: Re: Beware of NSolve - nastier example
- To: mathgroup at smc.vnet.net
- Subject: [mg50384] Re: [mg50372] Re: [mg50346] Re: Beware of NSolve - nastier example
- From: DrBob <drbob at bigfoot.com>
- Date: Wed, 1 Sep 2004 01:49:23 -0400 (EDT)
- References: <200408200858.EAA12533@smc.vnet.net> <cg6srb$odf$1@smc.vnet.net> <200408280837.EAA19074@smc.vnet.net> <200408311028.GAA18653@smc.vnet.net>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
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