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

```

• Prev by Date: RE: Gradient in FindMininum
• Next by Date: Re: Gradient in FindMininum
• Previous by thread: Bessel function / MeijerG / Integral bug?
• Next by thread: Re: Precision example