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

```

• Prev by Date: Re: Question: weight matrices and MonomialOrder for GroebnerBasis
• Previous by thread: Re: Beware of NSolve - nastier example
• Next by thread: 3D Screensavers for Mathematica users