Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: A problem of numerical precision

• To: mathgroup at smc.vnet.net
• Subject: [mg52642] Re: [mg52610] A problem of numerical precision
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Sun, 5 Dec 2004 02:08:38 -0500 (EST)
• References: <200412030815.DAA25377@smc.vnet.net> <200412040907.EAA13400@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Guofeng Zhang wrote:
>  
>  Hi,
> 
>  I met one problem when I do some iteration: By increasing numerical
>  precision, the results are so different from the original one!  I
>  don't know what went wrong, and hope to get some answers.
>  
>  The code is
>  
>  delta = 1/100;
>  a = 9/10;
>  b = -3*1.4142135623730950/10;
>  A = { {1,0}, {b,a} };
>  B= { {-1,1}, {-b,b} };
>  
>  f[v_,x_] := If[ Abs[ (10^20)*v-(10^20)*x]>(10^20)*delta, 1, 0 ];
>  
>  M = 3000;
>  it = Table[0, {i,M}, {j,2} ];
>  it[ [1] ] = { -delta/2, (a+b)*(-delta/2) };
>  
>  For[ i=1, i<M, it[ [i+1]]  = A.it[[i]]+f[ it[ [i,1] ], it[ [i,2] ]
>  ]*B.it[ [i] ]; i++ ];
>  temp = Take[it, -1000];
>  
>  ListPlot[ temp ];
>  
>  By evaluating this, I got an oscillating orbit. I got the same using
>  another system. However, if I increase the numerical precision by using
>  b = -3*1.4142135623730950``200/10;
>  to substitute the original b, the trajectory would converge to a fixed
>  point instead of  wandering around!
>  
>  I don't know why. I am hoping for suggestions. Thanks a lot.
>  
>  Guofeng
>  


You are seeing the effect of plotting "zeros" of negative accuracy.

Try it as below.

delta = 1/100;
a = 9/10;
b = -3/10*Sqrt[2];
A = {{1,0}, {b,a}};
B = {{-1,1}, {-b,b}};
f2[v_,x_] := UnitStep[Abs[v-x]-delta]
M = 3000;

itm = Table[0, {i,M}, {j,2}];
itm[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, MachinePrecision];
Do[itm[[i+1]] = A.itm[[i]] + f2[itm[[i,1]],itm[[i,2]]]*B.itm[[i]],
   {i,1,M-1}]

it200 = Table[0, {i,M}, {j,2}];
it200[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, 200];
Do[it200[[i+1]] = A.it200[[i]] + f2[it200[[i,1]],it200[[i,2]]]*B.it200[[i]],
   {i,1,M-1}]

it500 = Table[0, {i,M}, {j,2}];
it500[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, 500];
Do[it500[[i+1]] = A.it500[[i]] + f2[it500[[i,1]],it500[[i,2]]]*B.it500[[i]],
   {i,1,M-1}]

(These could also be computed with FoldList).

Now look at the last elements, or take the Precision, of it500. You will 
find that fewer than 13 digits are correct. At lower initial precision 
you get what amount to indeterminate values (except we have some 
information bounding their magnitudes). Using machine arithmetic there 
is no way to tell from the result whether it is "correct" or just how 
far off it might be. It turns out that they are just fine in comparison 
to it500. So plots of itm and it500 are reliable.

The last 1000+ values for it200 are a different matter. They all appear 
to be zero. For example:

In[156]:= it200[[1888]] == {0,0}
Out[156]= True

A closer inspection reveals that they have negative accuracy. The 
1888th, for example, could be anywhere up to around 262 digits in 
magnitude (and it gets worse as we move further along in it200).

In[157]:= InputForm[it200[[1888]]]
Out[157]//InputForm= {0``-262.03379085102995, 0``-261.84220176123233}

So they plot as zeros, but that's not a very reliable indication of the 
actual values.


Daniel Lichtblau
Wolfram Research

• Follow-Ups:
  • Re: Re: A problem of numerical precision
    • From: DrBob <drbob@bigfoot.com>

• References:
  • A problem of numerical precision
    • From: Guofeng Zhang <guofengzhang@gmail.com>