Re: Coaxing N[] to work
- To: mathgroup at smc.vnet.net
 - Subject: [mg73394] Re: Coaxing N[] to work
 - From: Scott Hemphill <hemphill at hemphills.net>
 - Date: Thu, 15 Feb 2007 05:01:40 -0500 (EST)
 - References: <equo6r$in0$1@smc.vnet.net>
 - Reply-to: hemphill at alumni.caltech.edu
 
<p at dirac.org> writes:
> Sometimes N[,] doesn't appear to work.  Like here:
> 
> 
> x = {2.0, 3.0, 5.0};
> A = { {6.0, 2.0, 1.0}, {2.0, 3.0, 1.0}, {1.0, 1.0, 1.0} };
> For[ k=0, k<15, ++k,
>    lambda = x.A.x/(x.x);
>    y = LinearSolve[A,x];
>    x = y / Norm[y,Infinity];
> ]
> N[lambda, 30]
> 
> 
> The output is:
> 
>    Out[5]= 0.578933
> 
> I was expecting 30 digits.  Why did N[] ignore my request for 30 digits?
Because lambda had already been calculated, and it didn't have 30 digits.
See how Mathematica prints the following numbers:
  In[1]:= 2
  Out[1]= 2
  In[2]:= 2.0
  Out[2]= 2.
  In[3]:= 2`30
  Out[3]= 2.00000000000000000000000000000
Now see what the precision of each number is:
  In[4]:= Precision[2]
  Out[4]= Infinity
  In[5]:= Precision[2.0]
  Out[5]= MachinePrecision
  In[6]:= Precision[2`30]
  Out[6]= 30.
Here's your problem, redone with infinite precision:
  In[7]:= x = {2, 3, 5};
  In[8]:= A = { {6, 2, 1}, {2, 3, 1}, {1, 1, 1} };
  In[9]:= For[ k=0, k<15, ++k,
             lambda = x.A.x/(x.x);
             y = LinearSolve[A,x];
             x = y / Norm[y,Infinity];
          ]
  In[10]:= lambda
           1102158619423970036337
  Out[10]= ----------------------
           1903774504398915184457
  In[11]:= N[lambda, 30]
  Out[11]= 0.578933385691052787623495851172
Here's your problem, redone with specified precision.  Whenever a number
with infinite precision is arithmetically combined with a number with
finite precision, the result has finite precision.  If I had wanted, I
could have set all the elements in x and A to finite precision by using
SetPrecision, e.g., "x = SetPrecision[ {2, 3, 5}, 30 ];".
  In[12]:= x = {2`30, 3, 5};
  In[13]:= A = { {6, 2, 1}, {2, 3, 1}, {1, 1, 1} };
  In[14]:= For[ k=0, k<15, ++k,
              lambda = x.A.x/(x.x);
              y = LinearSolve[A,x];
              x = y / Norm[y,Infinity];
           ]
  In[15]:= lambda
  Out[15]= 0.5789333856910527876234959
  In[16]:= Precision[%]
  Out[16]= 24.9696
Note that some precision was lost due to round off error.  You can't
get more precision after the number has already been calculated:
  In[17]:= N[%15,100]
  Out[17]= 0.5789333856910527876234959
Scott
-- 
Scott Hemphill	hemphill at alumni.caltech.edu
"This isn't flying.  This is falling, with style."  -- Buzz Lightyear