Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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