5
$\begingroup$

I need to calculate values of a highly nonlinear recursive function, and I am confused by the results Mathematica is returning.

z[n_, c_] := If[n > -1, z[n - 1, c]^2 + c, 0];
n = 11;
x = -Sqrt[2];
a1 = z[n, x]

$\left(\left(\left(\left(\left(\left(\left(\left(\left(\left(2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}\right)^2-\sqrt{2}$

a1//N

(* -0.046964 *)

This returned a value, but I am concerned about rounding errors in the innermost radical propagating through iterations and increasing. So I tried using Expand in square away some of the radicals, which worked:

a2=Expand[a1]

$147697092815735181686004274312885849093467177887186260181880117956815516746082505563737982134388389030893241923971655031156923982009795445237436098869911125858111671457534293531550227204426557474412874444056537558952359812249193088204461544748223088308707661946397504442822527292765990916280270193850801812897311597874449846714549460810335816823643227293330493863411548106839219241822466768921583568955627335268963141029610119304561256815146929152942585056781124212373816239809008225462890393635942384958615924327266051882051312940199368467608830911644642855433909402967926-104437615891545257496417768499303993150192546198507781167667629123050986331611062720622726045412096590344564060208770686303132757412581287221321140705193090702301237758123591562160043450173732009719642370018812330693084500590805964215749598013705892215307391548522404846842401369898675081900550147221730220216956830033181121514702889893933979494308972137475451707237696361816170276267908811592670531759321517526230879736820859410888074310814279229466095712432316237303167672249609959845163632221432738068908191253224393164837143730654199277516113353838271576786201217435197 \sqrt{2}$

Despite the large integers this form seems more usable because there's only one Sqrt[2] in it so a numerical value can be calculated with arbitrary and knowable precision. But,

a2//N

returns a "No significant digits are available to display" error.

So, what should I do here? Can a2//N be forced to calculate with sufficient precision? I don't know how a1 can be used without causing a snowballing error problem.

What I ultimately need is to plot n vs z[n,x] where n>200, so please check if proposed solutions work for large n. I use n=11 above because that's the lowest n at which I encounter this problem.

A related Question was asked 5 years ago, but the answered discuss Mathematica 8, I didn't entirely understand the accepted answer, and it doesn't directly address my issue of which method to use.

PS: If anyone knows how to fix my wonkily displayed outputs above, please do.

$\endgroup$
2
  • 1
    $\begingroup$ You can try Block[{$MaxExtraPrecision = 1000}, N[a2, 10] ]. $\endgroup$ Jan 10, 2016 at 11:35
  • $\begingroup$ Also SetPrecision[a2, 600], but you can get the same result with SetPrecision[a1, 30] as well. $\endgroup$
    – Mr.Wizard
    Jan 10, 2016 at 13:13

1 Answer 1

8
$\begingroup$

The loss of precision of this function by itself seems to be fairly modest: the big problem comes from the huge integers produced when expanding it. You can see that more clearly by changing the function to be tail recursive (so that larger values of n will be accessible without blowing the stack):

ClearAll[z2];
z2[result_: 0, n_Integer, c_] := z2[result^2 + c, n - 1, c];
z2[result_, -1, c_] := result;
SetAttributes[z, NHoldAll];

200 - Precision@z2[700, -Sqrt[2.`200]] (* -> 32.3568 lost digits *)

But actually this is overly conservative. And just as well, because N starts to choke when the expression gets very big (n about 700), taking exponentially longer for each increase in n...

200 + Log[10, N[z2[700, -Sqrt[2]], 200] - SetPrecision[z2[700, -Sqrt[2.`200]], 200]]
(* -> 13.4842 actual imprecise digits in the result *)

You'll have to increase $IterationLimit for n greater than 4093, but this is harmless. At that point, Mathematica thinks about 188 digits have been lost (although, as we know, they haven't really). I think it's probably easier and more reassuring to increase the precision of c than to SetPrecision the output of z2 before plotting.

As you can see, though, the difference in the actual values is fairly modest even for large n with machine-precision c:

Block[{$IterationLimit = Infinity},
  DiscretePlot[{z2[n, -Sqrt[2.`200]], z2[n, -Sqrt[2.`]]}, {n, 4050, 4150}]
 ]

DiscretePlot of 200-digit vs. machine-precision z2 output

The near-periodic behavior of z2[n, -Sqrt[2]] with respect to n for a several values near the $IterationLimit barrier seems interesting.

$\endgroup$
4
  • $\begingroup$ Hey, thank very much for this. I didn't understand some little pieces though. Could you explain what the ":" in your 2nd line means? (I don't know how to Google that). Also, just to be sure I understand, did you mean for the last argument of the 3rd line to be 0 instead of c? $\endgroup$ Jan 10, 2016 at 23:01
  • $\begingroup$ @JerryGuern the colon is Optional. The third argument on the third line is in the place where c would be and so takes the value c, but need not have any particular name as it's not used on the RHS anyway. I should correct one thing in the above: the tail-recursive form is not faster than the ordinarily recursive one. In fact, they are exactly the same. The advantage of tail recursion is just that you won't blow the stack when increasing n. $\endgroup$ Jan 10, 2016 at 23:41
  • $\begingroup$ Near-periodic it is. Too wobbly for any smallish periods. When I find values for which the period might be 4, or 8, it seems they are not attractive (in the sense that being epsilon away at a given step implies more than epsilon away 4 or 8 steps later, in fact around 3.4*epsilon). None of this is too surprising (at least after checking a couple of references); we are just inside the chaotic region for f(n)=f(n-1)^2+c. $\endgroup$ Jan 13, 2016 at 17:07
  • $\begingroup$ @DanielLichtblau Exactly. That's exactly what I'm beating MMa (and my MMa skills) up against, to make sure I know what I'm doing. $\endgroup$ Jan 16, 2016 at 22:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.