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)
- Beware of NSolve - nastier example