MathGroup Archive 2000

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

Search the Archive

Re: Precision example

  • To: mathgroup at smc.vnet.net
  • Subject: [mg24057] Re: Precision example
  • From: Mark Sofroniou <marks at wolfram.com>
  • Date: Thu, 22 Jun 2000 01:01:49 -0400 (EDT)
  • Organization: Wolfram Research Inc
  • References: <8ipn5k$an2@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

The reason for the discrepancy is that you are computing trigonometric
functions with very large arguments. The argument of Sin here grows
exponentially with the number of iterations:

Sin[2^n*ArcSin[Sqrt[x[0]]]]^2

When you compute Sin using double precision numbers, it simply
calls the C library math function sin for efficiency. Unfortunately,
some machine precision functions on some machines simply
return the argument unevaluated when it is very large - even if
that means that it is outside the appropriate range of the function.

This calls the C library machine double precision function sin.

In[1]:= Sin[8.181128596022229`*^35]

                  35
Out[1]= 8.18113 10

When you do the same computation in fixed precision software
arithmetic, then the bigfloat Sin function notices that there is
no possibility of obtaining a reasonable result from the given
input and so it returns 0.

In[2]:= prec = $MinPrecision = $MaxPrecision = 53/Log[2.,10.];

In[3]:= Sin[SetPrecision[8.181128596022229`*^35, prec]]

Out[3]= 0

Therefore after a certain point your exactSol returns 0 and
since your function x is positive, the difference using fixed
but arbitrary precision is always negative.

In our development version we have improved our build system to
detect platforms which return infeasible machine precision
results to avoid situations such as Out[1].

Mark Sofroniou,
Wolfram Research.



Arturas Acus wrote:

> 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: Conversion of Orderless functions to non Orderless one
  • Next by Date: Re: imposing side conditions on Solve
  • Previous by thread: Precision example
  • Next by thread: parallel computing toolkit