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