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