A difficult precision problem

Question

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.

Answer 1

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}]
 ]

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