Numerical instability of an extended tetration

Question

For bases $a\in(1,e^{1/e})$, ${}^na=a^{({}^{n-1}a)}=a^{a^{a^{.^{.^{.^a}}}}}$ converges to a value denoted as ${}^\infty a$. By observing the convergence rate of this sequence, we can derive the limit:

$$\lim_{n\to\infty}\frac{{}^\infty a-{}^{n+x}a}{{}^\infty a-{}^na}=[\ln({}^\infty a)]^x$$

By supposing we seek a continuous version of tetration that satisfies this, and rearranging so that ${}^xa$ is solved for, we derive:

$${}^xa=\lim_{n\to\infty}\log_a^{\circ n}({}^\infty a-({}^\infty a-{}^na)[\ln({}^\infty a)]^x)\tag1$$

where $\log^{\circ n}$ is the logarithm applied $n$ times. As an example, with $n=10$, I obtained the following plot:

which looks really nice. Then looking at $n=15$, I get this:

which raises concern. For $a$ close to $1$ and $n=10$, I get

It would seem to work well for small $n$ and large $a$, but then for larger $n$ or smaller $a$, it becomes unstable. As far as I can tell, this issue is due to the amount of numerical precision required while evaluating $(1)$, especially when the base is closer to $1$.

---

So the first question is whether this is due to numerical precision, or if it's simply because $(1)$ does not converge.

If it's the former, then is there any way to circumvent this without brute forcing with more precision? And how should I pick the values of $n$ for a given base $a$ (and $x$)?

If it's the latter, then does it converge anywhere?

---

Code for computing $(1)$, showing the following for $a=\sqrt2$ and $x=1.5$:

n   1.4142135623730^^n
--------------------------
0   1.42291711861386
1   1.4657586018199498
2   1.4910645646490854
3   1.5069501895748705
4   1.5172760309843982
5   1.5241342747726574
6   1.528753204049527
7   1.5318927292918296
8   1.5340399138955585
9   1.5355145848360043
10  1.5365302824374432
...
45  1.538805432574356
46  1.5388054445894592
47  1.5388054519338499
48  1.5388054652284342
49  1.5388054823911386
50  1.538805506512146
...
90  1.7233534923554696
91  1.755592017472159
92  2.0000000000000004
93  2.000000000000001
94  2.000000000000001

Showing the apparent value of $^{1.5}\sqrt2\simeq1.5388$ followed by divergence.

Answer 1

Conclusion of what follows below is that the formula/implementation is subject to catastrophic cancellation in the part ${}^\infty a-{}^na$.

I have implemented your code verbatim in MATLAB. The result with double precision is identical (subject to $\pm 1$ shift somehow in the counter $k$) for $k=25,50,100$.

>> test_double(25)

inf_a - n_a = 4.5947e-05

res =

   1.538798995072006


>> test_double(50)

inf_a - n_a = 4.8182e-09

res =

   1.538805596316793

>> test_double(100)

inf_a - n_a = 1.3323e-15  % Close to machine epsilon for 64 bits (!!!)

res =

   2.000000000000001

Now increasing the precision to quadruple (digits = 32)

>> test_vpa(25,32)

inf_a - n_a = 0.00004594710217926553685144416637951


res =

    1.538798995052591060169673385115

>> test_vpa(50,32)

inf_a - n_a = 0.0000000048181840417120357729198170741061


res =

    1.5388054167659345033199452184372

>> test_vpa(100,32)

inf_a - n_a = 0.00000000000000005298542667702117495721500428964


res =

    1.5388054174394144484626385875352

One can even go for more $k$, but for $k=200$ it becomes complex (maybe time to improve the bisection(?))

>> test_vpa(150,32)

inf_a - n_a = 0.00000000000000000000000058267893924374831695038060802142


res =

    1.5388052522455767035302332127693

>> test_vpa(200,32)

inf_a - n_a = -0.0000000000000000000000000000002180085870192506738614539008334


res =

    4.5484210435493524026130653949673 + 0.27008458814289299475281495540935i

---

Just for fun: $k=300$, 50 digits precision, bisection has 200 iterations:

>> test_vpa(300,50)

inf_a - n_a = 

ans =

0.00000000000000000000000000000000000000000000000077490532526797719918187015369561788340262562842819


res =

    1.5388054183673573750604847502205024901223348381524

Answer 2

Another example suggesting the problem is numerical instability of the double-precision based computations.

Using Pari/GP with internal precision of 200 decimal digits I get the following list:

    n  a^^1.5 (based on n iterations)
   ---------------------------------------
    1  1.4657586018199495028
    2  1.4910645646490847958
    3  1.5069501895748684802
    4  1.5172760309843957556
    .....
   27  1.5388009658025622487
   28  1.5388023318106759423
   29  1.5388032786497405906
   30  1.5388039349458743480
   ...
   83  1.5388054174394108432
   84  1.5388054174394125092
   85  1.5388054174394136640
   86  1.5388054174394144644
   87  1.5388054174394150192
   88  1.5388054174394154038
   89  1.5388054174394156704
   90  1.5388054174394158551
   91  1.5388054174394159832
   92  1.5388054174394160720
   93  1.5388054174394161335
   94  1.5388054174394161762
   95  1.5388054174394162057
   96  1.5388054174394162262
   97  1.5388054174394162404
   98  1.5388054174394162503
   99  1.5388054174394162571
  100  1.5388054174394162618
    ....
  200  1.538805417439416272500504006016378087864
  201  1.538805417439416272500504006016378484370
  202  1.538805417439416272500504006016378759207
  203  1.538805417439416272500504006016378949709
  204  1.538805417439416272500504006016379081755
  205  1.538805417439416272500504006016379173283
  206  1.538805417439416272500504006016379236725
  207  1.538805417439416272500504006016379280699
  208  1.538805417439416272500504006016379311180
  209  1.538805417439416272500504006016379332308

The Pari/GP-routine was

default(realprecision,200) \\ this is my standard internal precision request
\\ performs tetration to (i)nteger "heights" by base b
itet(x,h)=for(k=1,h,x=b^x);for(k=1,-h,x=log(x)/log(b));return(x)

b=sqrt(2)
\\ one could find the fixpoint (or "b^^infinity") by simple binary search
b_inf = solve(t=1,exp(1), b^t - t )
\\ but more efficient and more general is the following formula:
t0=b_inf=exp(-LambertW(-log(b)))  \\ fixpoint: compute b^^infinity

\\ tetration of an initial value z0 with base b
\\           to fractional heights x
\\               with n iterations, 
\\                      where n should go to infinity
{tetx(z0,x,max_n=10)=my(w);
   w= t0-(t0-itet(z0,max_n))*log(t0)^x;
   itet(w,-max_n)}

\\ compute a list of estimates based on n iterations, n=1..100
list = Mat(vectorv(100,n, [n, tetx(1,1.5,n)]))

\\ compute additional approximations, n=200.. 209
list1 = Mat(vectorv(10,n, [199+n, tetx(1, 1.5, 199+n )]))

Answer 3

We compare consecutive terms. Let $q=\ln({}^\infty a)$. Note that:

\begin{align}\log_a({}^\infty a-({}^\infty a-{}^{n+1}a)q^x)&=\log_a({}^\infty a)+\log_a\left(\frac{{}^\infty a-{}^{n+1}a}{{}^\infty a}q^x\right)\\&={}^\infty a+\frac{{}^\infty a-{}^{n+1}a}{\ln({}^\infty a)}q^x+\mathcal O(q^{2n})\\&={}^\infty a+({}^\infty a-{}^na)q^x+\mathcal O(q^{2n})\end{align}

Now we need to log both sides $n$ times. Note that

$$\log_a(\alpha+\epsilon)=\log_a(\alpha)+\frac{\epsilon}{\alpha\ln(a)}+\mathcal O(\epsilon^2)$$

We can also verify through a similar process that the limit is monotone increasing, and hence $\alpha\ge\alpha'$, which implies the error between consecutive terms to be $\mathcal O(q^2/\alpha'\ln(a))^n\simeq\mathcal O(q({}^\infty a/\alpha'))^n$. It is known that $q<1$, so all that remains to be seen is that $\alpha'$ is sufficiently close to ${}^\infty a$. From this, a rough estimate of $\alpha'$ can be deduced as ${}^xa$, with the error given from above, and hence for all $a$ and sufficiently large $x$, the limit exists, and furthermore, converges monotonically and uniformly to an analytic function.

By disregarding the smallest $N$ values that $\alpha$ takes, we can get a better bound of $\alpha'\simeq{}^{x+N}a$, and hence we can get the error's ratio to go below $1$ eventually, and hence this holds everywhere.

---

The above works for any $x\in\mathbb C$ since $[\ln({}^\infty a)]^{x+N}\to0$ as $N\to\infty$ even for complex $x$.