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

Precision example

• To: mathgroup at smc.vnet.net
• Subject: [mg24037] Precision example
• From: "Arturas Acus" <acus at itpa.lt>
• Date: Wed, 21 Jun 2000 02:19:59 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

Dear Group,
Recently there have been few letters on Precision and Accuracy of 
numerical calculations in Mathematica. 
The following example I think is worth to be explored.
Below is a recurency equation and it's exact solution. 
Unfortunately, not only Mathematica RSolve can't fint the exact 
solution, it but also such a 
clever algorithm as Hyper described in free internet book 
(http://www.cis.upenn.edu/wilf) 
was unable to find this particular solution.
In[2]:=
$Version
Out[2]=
"4.0 for Linux (April 21, 1999)"

One can verify this for example using FunctionExpand and Expand 
functions, for any integer n. Here is verification for n=3.
In[1]:=
x[n_Integer?(Positive[#] &)] := x[n] = 4*x[n - 1]*(1 - x[n - 1])
In[2]:=
Expand[x[3]]
Out[2]=
\!\(64\ x[0] - 1344\ x[0]\^2 + 10752\ x[0]\^3 - 42240\ x[0]\^4 + 
    90112\ x[0]\^5 - 106496\ x[0]\^6 + 65536\ x[0]\^7 - 16384\ 
x[0]\^8\)
In[3]:=
Expand[FunctionExpand[Sin[2^n*ArcSin[Sqrt[x[0]]]]^2 //. n -> 3]]
Out[3]=
\!\(64\ x[0] - 1344\ x[0]\^2 + 10752\ x[0]\^3 - 42240\ x[0]\^4 + 
    90112\ x[0]\^5 - 106496\ x[0]\^6 + 65536\ x[0]\^7 - 16384\ 
x[0]\^8\)
x[0] here is initial condition, chosen between -1 and 1 (due to Sin 
function). Let us make numerical simulations.
 Let us solve the recurrence equation numerically and then compare 
rezuls with
exact solution. First let's try machine arithmetic.

In[4]:=
x[0] = N[1/3];
x[n_Integer] := x[n] = 4*x[n - 1]*(1 - x[n - 1])
exactSol[n_Integer] := exactSol[n] = Sin[2^n*ArcSin[Sqrt[x[0]]]]^2
In[9]:=
ListPlot[Table[exactSol[n] - x[n], {n, 1, 70}], PlotRange -> {-1, 1}]
Out[9]=
Graphics

In[]:=Precision[exactSol[50] - x[50]]
Out[]:=16
So we can do safely only 48 iterations! After 48 iterations we see 
huge difference.
 Notice, that differences are positive and negative, 
thought I think there is more negative diferences.
Let us now quit kernel and make fixed precision float Mathematica  
simulation.
In[10]:=
Quit[]
In[1]:=
prec = $MinPrecision = $MaxPrecision = 35
Out[1]=
35
In[13]:=
x[0] = N[1/3, prec];
x[n_Integer] := x[n] = N[4*x[n - 1]*(1 - x[n - 1]), prec]
exactSol[n_Integer] := exactSol[n] = N[Sin[2^n*ArcSin[Sqrt[x[0]]]]^2, 
prec]
In[16]:=
ListPlot[Table[N[(exactSol[n] - x[n]), prec], {n, 1, 140}], 
  PlotRange -> {-1, 1}]
Out[16]=
Graphics
In[17]:=
Precision[exactSol[50] - x[50]]
Out[17]=
35

With increased fixed precision we can do safely 115 iterations 
roughtly. But what happened. We see only negative diferences in the plot!
  So Mathematica fixed precision
arithmetic is somehow different! Now let us do the same thing with 
Big Numbers arithmetic
In[18]:=
Quit[]
In[1]:=
prec = 35
Out[1]=
35
In[2]:=
x[0] = N[1/3, prec];
x[n_Integer] := x[n] = N[4*x[n - 1]*(1 - x[n - 1]), prec]
exactSol[n_Integer] := exactSol[n] = N[Sin[2^n*ArcSin[Sqrt[x[0]]]]^2, 
prec]
In[5]:=
ListPlot[Table[N[(exactSol[n] - x[n]), prec], {n, 1, 140}], 
  PlotRange -> {-1, 1}]
Out[5]=
Graphics

With Big Nums we can do only 60 iterations. This is not surprising 
due to precision loss in these simulations.
But now we can see that at 60-th iteration we have zero precision. 
Strange, but we see only positive diferences in Big Nums also.
In[6]:=
Precision [exactSol[60] - x[60]]
Out[6]=
0
In version 2.2 for Windows we could see both positive and negative 
diferences in these simulations. 
So something was changed in the behaviour since ver 2.2. Also it 
would
be interesting to know if these simulations are identical on other 
operating systems, say Windows. 
Of course, because Mathematica performs these simulations internaly, 
one could hope it should be so. 
Probably somebody from WRI can comment this issue, in particulary why 
in fixed precision Mathematica arithmetic 
simulation we see only negative differences, whereas true machine 
number simulations gives both?
                                     
Dr. Arturas Acus
Institute of Theoretical
Physics and Astronomy
Gostauto 12, 2600,Vilnius
Lithuania 


E-mail: acus at itpa.lt
   Fax: 370-2-225361
   Tel: 370-2-612906