Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Bug in NSolve?

• To: mathgroup at smc.vnet.net
• Subject: [mg33843] Re: [mg33788] Bug in NSolve?
• From: "Fred Simons" <f.h.simons at tue.nl>
• Date: Fri, 19 Apr 2002 02:29:46 -0400 (EDT)
• References: <200204160750.DAA20077@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

It looks like a bug, maybe related to some other strange behaviour I will
show in the examples. Maybe the problem has to do with small but not very
small imaginary parts of the solutions

We define a polynomial of degree depneding on two parameters such that the
imaginary parts of the two roots are-1/b and 1/b.

In[1]:=
Clear[a, b];
pol[a_, b_] = 1 + a^2 - 2*a*b*x + b^2*x^2;
Solve[pol[a, b] == 0, x]

Out[3]=
{{x -> (-I + a)/b}, {x -> (I + a)/b}}

In[4]:=
a = 45654; b = 1;
NSolve[pol[a, b] == 0, x]
pol[a, b] /. %

Out[5]=
{{x -> 45654. - 0.999999957203608*I},
  {x -> 45654. + 0.999999957203608*I}}

Out[6]=
{0. + 0.*I, 0. + 0.*I}

This is OK. But with one more digit in a we get incorrect results:

In[7]:=
a = 456544; b = 1;
NSolve[pol[a, b] == 0, x]
pol[a, b] /. %

Out[8]=
{{x -> 456543.9999666127}, {x -> 456544.0000333873}}

Out[9]=
{1.000030517578125, 0.999969482421875}

Fortunately, using arbitrary precision numbers helps:

In[10]:=
a = 456544; b = 1;
NSolve[pol[a, b] == 0, x, 20]
pol[a, b] /. %

Out[11]=
{{x -> 456543.999999999999999999771731`14.0141 -
     0.9999999999999999958`8.3546*I},
  {x -> 456543.999999999999999999771731`14.0141 +
     0.9999999999999999958`8.3546*I}}

Out[12]=
{-1.14844`0*^-9 + 9.2294`0*^-10*I,
  -1.14844`0*^-9 + -9.2294`0*^-10*I}

Now we look at b:

In[13]:=
a = 456544; b = 276355;
NSolve[pol[a, b] == 0, x, 20]
pol[a, b] /. %

Out[14]=
{{x -> 1.6520200466790903005192583234`14.014 -
     3.618534131823198421`8.3545*^-6*I},
  {x -> 1.6520200466790903005192583234`14.014 +
     3.618534131823198421`8.3545*^-6*I}}

Out[15]=
{-1.271611`0*^-8 + 9.90684`0*^-9*I,
  -1.271611`0*^-8 + -9.90684`0*^-9*I}

One digit more in b gives the strange and incorrect result observed by
rodgerro at library2.airnews.net :

In[16]:=
a = 456544; b = 2763556;
NSolve[pol[a, b] == 0, x, 20]
pol[a, b] /. %

Out[17]=
{{x -> 0.1652016459952322297793120673`14.0137}}

Out[18]=
{1.`2.0927}

One more increasing precision helps:

In[19]:=
a = 456544; b = 2763556;
N[NSolve[pol[a, b] == 0, x, 39]]
pol[a, b] /. %

Out[20]=
{{x -> 0.16520164599523224 - 3.6185262755666976*^-7*
      I}, {x -> 0.16520164599523224 +
     3.6185262755666976*^-7*I}}

Out[21]=
{0. + 0.*I, 0. + 0.*I}

One digit more in b yields incorrect results:

In[22]:=
a = 456544; b = 27635567;
N[NSolve[pol[a, b] == 0, x, 39]]
pol[a, b] /. %

Out[23]=
{{x -> 0.016520160415018807}}

Out[24]=
{0.99993896484375}

And so on?


Fred Simons
Eindhoven University of Technology

• References:
  • Bug in NSolve?
    • From: rodgerroSPAMNOT.siteconnect.com@library2.airnews.net