$\lim_{n\to\infty} f(2^n)$ for some very slowly increasing function $f(n)$

I should be able to answer this myself, but feel insecure anyway. I want to know, whether a function f(n) is bounded if n goes to infinity (and if it's bounded, the limit). Heuristically it appears (and I assume for generality) that f(n) is monotonuously increasing. But the increase of f(n) for n=1,2,3,4,... diminuishes and -if it eventually converges, then the convergence is much too slow for the estimate for some limit.
So I look at the sequence of f(2^n) and find, that the increase (="$\Delta(n)$") seems to halve at each step. Here is the table of values.

$\small \begin{array} {cccc} 2^n & f(2^n)& \Delta(n) =f(2^n)-f(2^n/2) & Q(n)=\Delta(n+1)/\Delta(n) \\ 2 & 1.00000000000 & 1.00000000000 & 1.00000000000 \\ 4 & 1.20710678119 & 0.207106781187 & 0.207106781187 \\ 8 & 1.32489552200 & 0.117788740818 & 0.568734351153 \\ 16 & 1.38378655210 & 0.0588910300981 & 0.499971641511 \\ 32 & 1.41323838275 & 0.0294518306501 & 0.500107242156 \\ 64 & 1.42796608046 & 0.0147276977050 & 0.500060518479 \\ 128 & 1.43533039941 & 0.00736431894937 & 0.500031919236 \\ 256 & 1.43901267941 & 0.00368228000046 & 0.500016366181 \\ 512 & 1.44085384991 & 0.00184117050297 & 0.500008283655 \\ 1024 & 1.44177444283 & 0.000920592923335 & 0.500004166833 \\ 2048 & 1.44223474122 & 0.000460298385385 & 0.500002089651 \\ 4096 & 1.44246489089 & 0.000230149674341 & 0.500001046382 \\ 8192 & 1.44257996585 & 0.000115074957672 & 0.500000523580 \end{array}$

As Q(n) tends to 0.5, I think I should conclude, that f(n) is bounded because the sum of all $\Delta$ is bounded, but I feel not sure, whether this is a meaningful conclusion here.
And if it were meaningful: what would then be a reasonable estimate for the upper bound?

[update]: Heuristics suggest, that $$f(n) \approx { {n-1 \over \log(2)} + {1 \over 2} \over n}$$ and then $$\lim_{n \to \infty} f(n) = {1 \over \log(2)} \approx 1.44269504089$$
[end update]

The function results of an adaption of the problem, whether for some integer exponent n the sum of the n'th power of the first m numbers can equal (m+1)^n, say whether there exist some natural n>2, m>2, such that $S(m,n) = (m+1)^n$ where $S(m,n)=\sum_{k=1}^m k^n$ .

I convert this in a formulation using the bernoulli-polynomials for the sum of the n'th powers. It is not difficult to find, that m must be greater than n.
By use of the bernoulli-polynomials one can also meaningfully interpolate m to fractional values, say x, in the sum S(m,n). So d(x,n) is the difference $d(x,n) = B(x,n)- (x+1)^n$ where B(x,n) is the Bernoulli-polynomial-generalized form of S(m,n).

With this, for a given integer n, I solve for x such that d(x,n)=0. Again, in general x>n but I want to know, whether x/n is bounded... So I determine $f(n)={\text{<root x of d(x,n)>} \over n}$
Remark: in my earlier notation above in the table I used x where I meant the exponent which I now more meaningfully denote as n.

Well, while I'm writing this: possibly the current problem can be solved on a much simpler way.... (At least I've now a meaningful heuristic, see above)

But my deeper curiosity is about the idea; whether -and to what extend- the above idea of "acceleration of convergence" is/can be made valid for such a type of conclusion at all?

-
 What are the assumptions on $f$? Better yet would be if you gave the function $f$. With just statistical data, its impossible to say anything with certainty. – picakhu Apr 30 '11 at 5:23 @picakhu: well I thought, the table would allow a general statement; I'll append a description of the function in the text of the question – Gottfried Helms Apr 30 '11 at 5:28 I suspect you could use something like Richardson extrapolation here (if $f(2^n)$ really does approach a limit). The error ratio $Q(n)$ seems to be converging to 0.5 so I would think this could be a candidate for extrapolation... – J. M. Apr 30 '11 at 7:57 Gottfried, I might have time to play with this later; could you share the code you used for crunching those numbers? – J. M. Apr 30 '11 at 9:41 @J.M.: Ok; I'll put the Pari/GP code into a new answer-field because the comment-boxes might be too small for some expression. I'll need a bit time... – Gottfried Helms Apr 30 '11 at 10:30

This is more of a (numerical) plausibility argument to support Gottfried's $\lim\limits_{n \to \infty} f(n) = \frac1{\log\;2}$ conjecture than an answer.

What struck me was that the error ratios $Q(n)$ as $n$ is doubled seem to indicate a decrease in error by a factor of 2. This reminded me of the decrease in error exhibited by the trapezoidal rule for numerical integration as the number of subintervals is doubled, leading me to think that Richardson extrapolation might be appropriate to use in numerically approximating the limit, just like the way Richardson extrapolation is used to improve trapezoidal estimates in Romberg quadrature.

Recall that what Richardson extrapolation essentially does is to fit an interpolating polynomial $p(x)$ to successive values of a sequence $s_i$, along with some auxiliary sequence $x_i$ that tends to 0 as $i\to\infty$ such that $p(x_i)=s_i$, and then take $p(0)$ to be the estimate of the limit of the sequence $s_i$. For this particular case, the error ratio appearing to approach 2 suggests that we take the auxiliary sequence $x_j=2^{-j}$.

The Richardson scheme proceeds as follows: letting $T_j^{(0)}=s_j$, we generate a triangular array of values

$$\begin{matrix}T_0^{(0)}&&&&\\T_1^{(0)}&T_1^{(1)}&&&\\T_2^{(0)}&T_2^{(1)}&T_2^{(2)}&&\\\vdots&\vdots&\vdots&\ddots&\\T_n^{(0)}&T_n^{(1)}&\cdots&\cdots&T_n^{(n)}\end{matrix}$$

whose entries are generated according to the recursive formula

$$T_j^{(n)}=T_{j}^{(n-1)}+\frac{T_{j}^{(n-1)}-T_{j-1}^{(n-1)}}{\frac{x_{j-n}}{x_j}-1}$$

and the $T_n^{(n)}$ (i.e., the "diagonal" elements of the array) form a sequence which (hopefully) converges faster to the limit than the original $s_j$. For this specific example, we have $s_j=f(2^{j+1})$ and $\frac{x_{j-n}}{x_j}=2^n$.

The Richardson scheme can be implemented such that only a one-dimensional array is needed, and this is what I will do in the Mathematica code that follows (specializing to the case $\frac{x_{j-n}}{x_j}=2^n$):

richardson[seq_?VectorQ] := Module[{n = Length[seq], m, res, ta},
ta[1] = First[seq]; res = {First[seq]};
Do[
ta[k + 1] = seq[[k + 1]];
m = 1;
Do[
m *= 2;
ta[j] = ta[j + 1] + (ta[j + 1] - ta[j])/(m - 1);
, {j, k, 1, -1}];
res = {res, ta[1]};
, {k, n - 1}];
Flatten[res]
]


We then generate the (first fourteen members of the) sequence in the OP to 100 significant figures:

pa[n_Integer, x_] := Expand[Sum[(-1)^k*BernoulliB[k]*Binomial[n, k]*
(x^(n - k + 1)/(n - k + 1)),
{k, n}] +  x^(n + 1)/(n + 1) - (x + 1)^n]

vals = Table[x/2^j /. FindRoot[pa[2^j - 1, x] == 0, {x, 2^j - 1, 2^(j + 1) - 2},
WorkingPrecision -> 100], {j, 14}];

Last[vals] (* sample value, j == 14 *)
1.442637503359895400534264282420018712540776326004454785114823591142583865448308469439211266125019043


If we apply richardson[] to vals, we get the following:

Last[richardson[vals]]
1.44269504088896340735992467998804478392151992973760702424401781700078935678750451268434910460102062

% - 1/Log[2]
-1.0138473535051260244153789098914315899303198623936805672011775182924857213751333727278493*^-27


and the result of the extrapolation matches $\frac1{\log\;2}$ to ~ 27 digits.

Generating more members of the sequence and computing with more digits might give a more accurate value, but instead of doing this, I decided to apply a different extrapolation algorithm, the Shanks transformation, with the reasoning that if two unrelated algorithms give (almost) the same answer for a given sequence, the answer is likely correct.

I will skip the details of the mathematics behind the Shanks transformation (but see this article) and instead present Mathematica code:

wynnEpsilon[seq_?VectorQ] := Module[{n = Length[seq], ep, res, v, w},
res = {};
Do[
ep[k] = seq[[k]];
w = 0;
Do[
v = w; w = ep[j];
ep[j] =
v + (If[Abs[ep[j + 1] - w] > 10^-(Precision[w]), ep[j + 1] - w,
10^-(Precision[w])])^-1;
, {j, k - 1, 1, -1}];
res = {res, ep[If[OddQ[k], 1, 2]]};
, {k, n}];
Flatten[res]
]


This is essentially the algorithm behind the hidden built-in Mathematica function SequenceLimit[], but I have chosen to implement the algorithm explicitly for pedagogical reasons.

We then have the results

Last[wynnEpsilon[vals]]
1.442695040888963407376933999717472189155804331509810458775727421372220122967729979363328849794520

% - 1/Log[2]
1.70093187155800517291583773568245246402780144411109037865448994778022269010133063583652*^-20


with a result that matches $\frac1{\log\;2}$ to ~ 20 digits. (The slight loss of accuracy is the price for the generality of the Shanks transformation, which did not exploit the behavior of the error ratios.)

The results of these two different extrapolation algorithms lends credibility to the conjecture that the limit is $\frac1{\log\;2}$.

-
Very nice, thank you. These examaples should help me to understand that extrapolation-procedures. – Gottfried Helms Apr 30 '11 at 21:16
The two extrapolation methods I mentioned in this post are discussed in Appendix A of the book The SIAM 100-Digit Challenge, and in more detail in Brezinski/Redivo Zaglia's Extrapolation Methods: theory and practice. – J. M. May 1 '11 at 0:58
I'll "accept" this answer to close the case. Many thanks to J.M. for the elaboration of the Richardson/Shanks-procedures which give me unexpected extra-profit...;-) <hr> Meanwhile I've seen the 1/log(2)-value in some of the literature where I was hinted to by the other thread. The question, whether that boundedness-criterion is valid seems now nearly trivial to me after we can insert a constant (but even bigger) value for the quotient to get a bounded sum for the deltas - what made me insure was the increasing stepwidth of the x - but this cannot be a problem. Thanks to all! – Gottfried Helms May 4 '11 at 7:48

Conceivably, the function $f(x)$ is within $0.00005$ of $1.44$ up until $x + 2^{2^{2^{10}}}$, but only then does the function suddenly skyrocket. Knowing localized data doesn't tell us anything about the limits of the function without some other sort of knowledge, such as the form of the function, the characteristics of its variation or derivative, etc.

-
 Not only that, but knowing something about the value of $f$ at powers of $2$ tells you nothing about what it does in between. Still, I think it might be helpful if someone could answer this: suppose $f$ is monotone increasing, and $Q(x)$ (as defined in the question) is monotone decreasing to $1/2$. Does that imply boundedness (hence, convergence) of $f$? – Gerry Myerson Apr 30 '11 at 5:46 Yes, f(x) surely is monotonously increasing. I added the description of f in my question – Gottfried Helms Apr 30 '11 at 6:11

This is not an answer, but Pari/GP-code for the computation of the data according to J.M.'s comment.
For floatingpoint computation I use internal precision of 200 dec digits per default.

\\ ==================== Pari/GP-Code =========================================
\\ Bernoulli-polynomials (in the Faulhaber-form, expressed in terms of zeta)
\\           keep exponent in global variable bn
{ init_bp(n) = bn = n;            \\ set global variable bn for exponent
bp = vector(2+n);       \\ coefficients of bernoulli-polynomial in global variable bp
bp[1+0] = 0;          \\ remember: Pari/GP indexes begin with offset 1
for(c = 1, n, bp[1+c] = -zeta(c-n) * binomial(n,c)) ;
bp[1+(1+n)] = 1/(1+n);
return(bp); }

\\ sum-of-like-powers 1^n + 2^n +... + x^n
\\     by bernoulli-polynomials
solp(x) = return( sum(k = 0, #bp -1, x^k * bp[1+k] )  )

\\ difference-function   (1^n+2^n+...+x^n) - (x+1)^n
solp_d(x)= solp(x) - (x+1)^bn

\\ final criterion-function : find root for some solp_d(n)
\\ solp_d(n)/(n+1)
{ f(n) = init_bp(n);
x0 = solve(x = bn, 2*bn, solp_d(x) ) ;   \\ find root between n and 2n
return( x0 /(bn+1) ); }

\\ -----------------------------------------------------
\\ remark: in the documented data of my question I had an additional
\\ index-offset of 1, so the data for displayed exponent n=8 were
\\ really for n=7
\\ so disp_f(512) gives the same value as documented for 512
disp_f(n) = f(n-1)

disp_f(128) \\ should be 1.43533039941

\\ -----------------------------------
delta(n)= if(n==1, disp_f(2^n), disp_f(2^n)-disp_f(2^n/2))

delta(2)  \\ should be 0.207106781187

\\ -----------------------------------
q(n) = if(n==1, 1, delta(n)/delta(n-1))

q(3) \\ should be 0.568734351153

-
 A question: solp(x) is of degree n+1? – J. M. Apr 30 '11 at 12:40 @J.M.: yes (..... dummytext to allow commentformat) – Gottfried Helms Apr 30 '11 at 12:51 Just to check my implementation: if n is 4, solp(x) is $\frac{x^5}{5}+\frac{x^4}{2}+\frac{x^3}{3}-\frac{x}{30}$? – J. M. Apr 30 '11 at 12:54 Yepp. this is the same what I have. – Gottfried Helms Apr 30 '11 at 12:58 After experimenting with both Richardson extrapolation and the Shanks transformation, the numerical results I'm getting seems to suggest that your heuristic of $\frac1{\log\;2}$ as the limit is accurate. If you're interested, I can post my experiment... – J. M. Apr 30 '11 at 14:05