Convergence acceleration for a special limit

Question

This question is inspired by another one on MSE: if a sequence $x_n$ is defined by $x_1 = \sqrt{2}, \quad x_{n+1} = x_n + \frac{1}{\sqrt{n} \cdot x_n} \quad \text{for } n \ge 1$, how to prove that it diverges?
Squaring the recurrence gives $$x^2_{n+1}=x^2_n+\frac2{\sqrt{n}}+\frac1{n x^2_n},\tag1$$ so we have
$$x^2_{n+1}>x^2_n+\frac2{\sqrt{n}}>x^2_n+\frac4{\sqrt{n}+\sqrt{n+1}}=x^2_n+4(\sqrt{n+1}-\sqrt{n}),$$ i.e. $$x^2_{n+1}-4\sqrt{n+1}>x^2_n-4\sqrt{n}\ge x^2_1-4=-2.$$ This means that $x^2_n-4\sqrt{n}$ is monotone increasing. From (1), we have $$x^2_{n+1}-4\sqrt{n+1}=(x^2_n-4\sqrt{n})+\frac2{\sqrt{n}(\sqrt{n}+\sqrt{n+1})^2}+\frac1{n x^2_n},$$ because $$\frac2{\sqrt{n}}-4(\sqrt{n+1}-\sqrt{n})=\frac2{\sqrt{n}(\sqrt{n}+\sqrt{n+1})^2}.\tag2$$ With $y_n=x^2_n-4\sqrt{n}$, we have $$y_{n+1}-y_n=\frac2{\sqrt{n}(\sqrt{n}+\sqrt{n+1})^2}+\frac1{n (4\sqrt{n}+y_n)}\tag3.$$ Since the RHS is $O(n^{-3/2})$, $$y_\infty-y_1=\lim_{n\to\infty}y_n-y_1=\sum^\infty_{n=1}\left(\frac2{\sqrt{n}(\sqrt{n}+\sqrt{n+1})^2}+\frac1{n (4\sqrt{n}+y_n)}\right)<\infty\tag4.$$

How can we calculate $y_\infty$ accurately? And: how can we 'accelerate convergence' for this limit?

From (4), it's clear that $y_\infty-y_n=O(n^{-1/2})$, just a few digits for quite some effort. Some acceleration of convergence is possible, though:
Using the known fact $$\lim_{n\to\infty}\left(\sum^n_{k=1}\frac1{\sqrt{k}}-2\sqrt{n}\right)=\zeta\left(\frac12\right)$$ and (2), we see that $$\sum^\infty_{n=1}\frac2{\sqrt{n}(\sqrt{n}+\sqrt{n+1})^2}=4+2\zeta\left(\frac12\right).$$ The second ingredient of the RHS of (4) doesn't converge any better, but we can extract the asymptotically dominant part: $$\frac1{n (4\sqrt{n}+y_n)}=\frac1{4n\sqrt{n}}-\frac{y_n}{4n\sqrt{n} (4\sqrt{n}+y_n)},$$ Thus, we obtain $$y_\infty-y_1=4+2\zeta\left(\frac12\right)+\frac14\zeta\left(\frac32\right)-\sum^\infty_{n=1}\frac{y_n}{4n\sqrt{n} (4\sqrt{n}+y_n)}.\tag5$$ Of course, we have to approximate the value with a finite sum: $$y_\infty\approx2+2\zeta\left(\frac12\right)+\frac14\zeta\left(\frac32\right)-\sum^{N-1}_{n=1}\frac{y_n}{4n\sqrt{n} (4\sqrt{n}+y_n)}.$$ Unfortunately, summing that for larger values of $N$ can't be done in double precision, because then, $1+10^{-16}=1$ after rounding. Working in higher precision gives $y_\infty\approx\underline{0.02230651}497489885274934240454627263466$ for $N=10^8$, the underlined digits are presumably correct.
In principle, we can continue that: $$\frac{y_n}{4n\sqrt{n} (4\sqrt{n}+y_n)}=\frac{y_n}{16n(n+1)}-\frac{y_n(\sqrt{n}y_n-4)}{16n(n+1)\sqrt{n}(4\sqrt{n}+y_n)}.$$ The second term on the RHS is $O(n^{-5/2})$, while summation by parts gives $$\sum^\infty_{n=1}\frac{y_n}{16n(n+1)}=\frac{y_1}{16}+\sum^\infty_{n=1}\frac{y_{n+1}-y_n}{16(n+1)}.$$ We can plug in the value of $y_{n+1}-y_n$ from (3), and using the first $10^8$ terms in each series gives $$y_\infty\approx\underline{0.0223065149610}5103102817029712459592728.$$ Now my question is: can we calculate the limit with significantly higher accuracy, but without still more hideous formulas?

Answer 1

Substituting the ansatz $$ y_n \sim \sum_{k = 0}^\infty \frac{c_k}{n^{k/2}} , \quad c_0=y_{\infty}, $$ into the recurrence relation and expanding in powers of $n^{-1/2}$ produces a recurrence for the coefficients $c_k$. The first few coefficients are \begin{align*} c_1 &= -\tfrac{3}{2}, \\[6pt] c_2 &= \tfrac{y_{\infty}}{16}, \\[6pt] c_3 &= \tfrac{-26 - y_{\infty}^2}{96}, \\[6pt] c_4 &= \tfrac{y_{\infty}\left(29 + y_{\infty}^2\right)}{512}, \\[6pt] c_5 &= \tfrac{-760 - 122y_{\infty}^2 - 3y_{\infty}^4}{7680}, \\[6pt] c_6 &= \tfrac{y_{\infty}\left(3655 + 317y_{\infty}^2 + 6y_{\infty}^4\right)}{73728}, \\[6pt] c_7 &= \tfrac{-14200 - 5618y_{\infty}^2 - 327y_{\infty}^4 - 5y_{\infty}^6}{286720}, \\[6pt] c_8 &= \tfrac{y_{\infty}\left(985655 + 162718y_{\infty}^2 + 7047y_{\infty}^4 + 90y_{\infty}^6\right)}{23592960}, \\[6pt] c_9 &= \tfrac{-11639840 - 8087956y_{\infty}^2 - 841686y_{\infty}^4 - 28818y_{\infty}^6 - 315y_{\infty}^8}{371589120}, \\[6pt] c_{10} &= \tfrac{y_{\infty}\left(450436865 + 125437246y_{\infty}^2 + 9385017y_{\infty}^4 + 264402y_{\infty}^6 + 2520y_{\infty}^8\right)}{13212057600}, \\[6pt] c_{11} &= \tfrac{-2872104320 - 2938510180y_{\infty}^2 - 487160862y_{\infty}^4 - 28191798y_{\infty}^6 - 672345y_{\infty}^8 - 5670y_{\infty}^{10}}{130799370240}. \end{align*} These coefficients can be used to approximate $y_n$ for large $n$, and numerical methods can then be applied to find an approximate value of $y_{\infty}$ by solving the resulting polynomial equation.

For example, using all terms up to and including $\frac{c_{11}}{n^{11/2}}$ with $n=10^8$, I obtained $$ y_\infty = 0.0223065149610197860946440297805177788\ldots $$ with all displayed digits being correct.

Addendum. Define $$ b_k^{(1)} = \begin{cases} 4\binom{1/2}{k/2+1}, & k \text{ even},\, k \ge 2, \\ 0, & \text{otherwise}, \end{cases} $$ $$ b_k^{(2)} = \sum_{j=1}^{\lfloor k/2 \rfloor} c_{k-2j+1} \binom{j-(k+1)/2}{j}, $$ and $s_k = b_k^{(1)} + b_k^{(2)}$. Then $c_0=y_\infty$, $c_1=-\frac{3}{2}$, and $$ 2k\,c_k= 4\left(b_{k+1}^{(1)} + \sum_{j=2}^{\lfloor(k+1)/2\rfloor} c_{k-2j+2} \binom{j-k/2-1}{j}\right) + \sum_{r=1}^{k} c_{r-1}s_{k-r+1} $$ for $k\ge 2$.

Wolfram Mathematica code:

y = Symbol["y"];

b1[k_] := If[EvenQ[k] && k >= 2, 4 Binomial[1/2, k/2 + 1], 0];
b2[k_] := Sum[c[k - 2 j + 1] Binomial[j - (k + 1)/2, j], {j, 1, Floor[k/2]}];
s[k_] := b1[k] + b2[k];
c[0] = y;
c[1] = -3/2;
Do[c[k] = (4 (b1[k + 1] +
Sum[c[k - 2 j + 2] Binomial[j - k/2 - 1, j], {j, 2,
Floor[(k + 1)/2]}]) +
Sum[c[r - 1] s[k - r + 1], {r, 1, k}])/(2 k);
c[k] = Simplify[c[k]], {k, 2, 11}];
Table[c[k], {k, 0, 11}]

Answer 2

Please consider this post as a long comment :)

To the nice posted solutions we can also add a heuristic approach allowing us to get the asymptotics of $x_n$

Starting from the equation $$x_{n+1}-x_n=\frac1{\sqrt n x_n}$$ choosing $x_n=2n^{1/4}+\epsilon_n$ and using the Taylor series $\epsilon_{n+1}-\epsilon_n=\frac{d\epsilon_n}{dn}+\frac12\frac{d^2\epsilon_n}{dn^2}+...$ $$(\epsilon_{n+1}-\epsilon_n)+2\Big((n+1)^{1/4}-n^{1/4}\Big)=\frac1{\sqrt n\,(2n^{1/4}+\epsilon_n)}$$ what gives $$\frac{d\epsilon_n}{dn}+\frac12\frac{d^2\epsilon_n}{dn^2}+...+2n^{1/4}\left(\frac1{4n}-\frac 3{32n^2}+...\right)=\frac1{2n^{3/4}}-\frac{\epsilon_n}{4n}+\frac{\epsilon_n^2}{8n^{5/4}}+...$$ Supposing that $\epsilon_n=\epsilon_1+\epsilon_2+\epsilon_3+...\,$, where $\epsilon_3\ll\epsilon_2\ll\epsilon_1$, and using the step-ny step approach, we get the recurrence relation for $\epsilon_k\,$ : $$\frac{d\epsilon_1}{dn}=-\frac{\epsilon_1}{4n}\,\, \Rightarrow\,\,\epsilon_1=\frac{C}{n^{1/4}}$$ where $C$ is the constant of integration, depending on boundery conditions (on $x_1$ in our case) $$\frac{d\epsilon_2}{dn}+\frac12\frac{d^2\epsilon_1}{dn^2}-\frac3{16n^{7/4}}+...=-\frac{\epsilon_2}{4n}+\frac{\epsilon_1^2}{8n^{5/4}}+...$$ from where we get the equation for $\epsilon_2\,$ : $$\frac{d\epsilon_2}{dn}+\frac{\epsilon_2}{4n}=\left(\frac{C^2}8+\frac3{16}\right)\frac1{n^{7/4}}\,\,\Rightarrow\,\,\epsilon_1=-\frac{2C^2+3}{8\,n^{3/4}}$$ Hence, $$x_n=2n^{1/4}+\frac C{n^{1/4}}-\frac{2C^2+3}{8\,n^{3/4}}+...=2n^{1/4}\left(1+\frac C{2\sqrt n}-\frac{2C^2+3}{16n}+...\right)$$ where $C=\frac{y_\infty}4$ in the notations above. We cannot find this constant in the frame of purely asymptotic approach, as we need the bridge to link the boundary conditions ($x_n$ at small $n$) to the solution valid at $n\gg1$.

Answer 3

Finally, I've implemented a variant of @Gary's method using a package for handling formal power series. It can be used to generate hundreds of terms of the asymptotic expansion automagically, without ugly recurrences for the coefficients.

Starting point is another consequence of equation (3) in the question: $$y_\infty-y_n=\sum^\infty_{i=n}\left(\frac2{\sqrt{i}(\sqrt{i}+\sqrt{i+1})^2}+\frac1{i (4\sqrt{i}+y_i)}\right)\tag6$$ If we have an asymptotic expansion of $y_n$ as a formal power series in powers of $n^{-1/2}$ (formal, because asymptotic series usually don't converge), the RHS could be expanded, too, and we'd have to compute expressions like $\displaystyle\sum^\infty_{i=n}i^{-m/2}$, where $m\ge3$. But that can be done (asymptotically!) via the Euler–Maclaurin formula: $$\sum^\infty_{i=n}i^{-m/2}=\frac2{m-2}n^{-m/2+1}+\frac12n^{-m/2}-\sum^\infty_{k=1}b_kc_{m,k}n^{-m/2-2k+1},$$ where $b_k=\displaystyle\frac{B_{2k}}{(2k)!}$ and $c_{m,k}=(-m/2)_{2k-1}$ is the falling factorial. If we understand all of that as formal power series in $x=n^{-1/2}$, this amounts to a linear series transformation $$S\left(\sum^\infty_{m=3}a_mx^m\right)=\sum^\infty_{m=1}s_mx^m,$$ where $$s_m=\frac2ma_{m+2}+\frac12a_m-\sum^{\lfloor(m-1)/4\rfloor}_{k=1}b_kc_{m-4k+2,k}a_{m-4k+2}.$$ That's essentially the only thing we have to code, and then, we can write (6)(assuming $y=y(x)=y(n^{-1/2})=y_n$) as $$y=y_\infty-S\left(\frac{2x^3}{(1+\sqrt{1+x^2})^2}+\frac{x^3}{4+xy}\right).\tag7$$ It's not a lot of new code: $S$ is implemented as a function tailsum:

fun tailSum(t: Taylor) = Taylor({ n -> if (n < 1) 0.0 else 2.0 / n * t.coeff(n + 2) + 0.5 * t.coeff(n) - (1..(n - 1) / 4).map { k -> b(k) * c(n - 4 * k + 2, k) * t.coeff(n - 4 * k + 2)}.sum() })

Yes, I know, functional programming, matter of taste.

The implementation of (7) is incredibly simple:

            val y = Taylor()
        y.assign(yInf - tailSum(shiftRight(2.0 / (2.0 + x * x + 2.0 * sqrt(1.0 + x * x)) + 1.0 / (4.0 + x * y), 3)))

So for any $y_\infty$, we get an asymptotic expansion, and we can compare the computed value of $y_n$ with the asymptotic value, using $2n$ terms of the series (nearly optimal):

d = y.value(2 * n, 1.0 / sqrt(n.toDouble())) - yN

Simple iteration would just mean

yInf -= d

and that converges quite rapidly, even when using $y_n$ as a starting value. In double precision, we get $y_\infty=0.02230651496101979$ in $7$ iterations (using $n=40$).

Can Newton iteration accelerate convergence? We don't calculate explicit formulas for the coefficients in terms of $y_\infty$, as Gary did, but that's not really necessary: with $\displaystyle z=\frac{\partial y}{\partial y_\infty}$, from (7) follows at once $$z=1+S\left(\frac{x^4z}{(4+xy)^2}\right),$$ so all we have to add to the code is

        z.assign(1.0 + tailSum(shiftRight(z / pow((4.0 + x * y), 2.0), 4)))

      
   yInf -= d / z.value(2 * n, 1.0 / sqrt(n.toDouble()))

Then, we get $y_\infty=0.02230651496101978$ in $4$ iterations.

For really higher precision, I had to adapt the power series code to use multiprecision coefficients, but otherwise, the algorithm didn't change. To cut a long story short: for $n=240$, using $480$ terms of the asymptotic expansion, we obtain $$y_\infty=0.0223065149610197860946440297805177788248978000325902772630432627686614713872729502390487156786893306$$ in 6 Newton steps, within a few seconds. Simple iteration needs 29 steps, roughly doubling the runtime, but giving the same result, naturally.

Strange fact: while this method gives consistent results with more and more digits while increasing $n$, they differ from Gary's result already after the 25th decimal. At first, I suspected a bug in my homespun power series/multiprecision code, so I made a check, computing $$\zeta(1/2)=\lim_{n\to\infty}\left(\sum^n_{i=1}\frac1{\sqrt{i}}-2\sqrt{n}\right)$$ with the same method, same $n$ and $2n$ terms. The result was the same as that calculated with a completely different method, and Pari/GP disagreed only concerning the 100th decimal. To resolve that, I plan to port Gary's code to Pari/GP and to compare results.

EDIT: the discrepancy is resolved, mea culpa. The value of $y_{10^8}$ I delivered in a comment was wrong, off by about $1.38e-26$, no idea how that could happen. Rounding errors are inevitable if you perform about $1e8$ operations, but that wouldn't explain an error of that magnitude. Checks with my own library or Pari/GP weren't much better, so I re-calculated it with a bigger working precision. Rounded to $38$ decimals, that gives $$y_{10^8}=0.02215651497469051955417471310030332440.$$ With that value, and Gary's asymptotics (and temporarily increasing working precision), that would give $$y_{\infty}=\underline{0.0223065149610197860946440297805177788}167483127,$$ which is compatible with my results.

Answer 4

For the divergence.

$$ x_{n+1} = x_n + \frac{1}{\sqrt{n} \, x_n} $$

In a first step, write $$ x_{n+1}- x_n =\frac{1}{\sqrt{n} \, x_n} $$ which is "close" to $$x'(n)=\frac{1}{\sqrt{n} \, x(n)}\quad \implies \quad x(n)= \sqrt{4 \sqrt{n}-2}$$ which is almost decent.

$$\left( \begin{array}{cccc} n & x_n & x(n) & x_n-x(n) \\ 2 & 2.12132 & 1.91229 & 0.209030 \\ 3 & 2.45465 & 2.21996 & 0.234698 \\ 4 & 2.68986 & 2.44949 & 0.240370 \\ 5 & 2.87574 & 2.63520 & 0.240545 \\ 6 & 3.03126 & 2.79248 & 0.238773 \\ 7 & 3.16594 & 2.92968 & 0.236259 \\ 8 & 3.28532 & 3.05184 & 0.233483 \\ 9 & 3.39294 & 3.16228 & 0.230659 \\ 10 & 3.49118 & 3.26330 & 0.227882 \\ 20 & 4.19195 & 3.98604 & 0.205904 \\ \end{array} \right)$$ and $$\left( \begin{array}{ccc} n & \frac {x_{n+1}}{x_n} & \frac {x(n+1)}{x(n)}\\ 2 & 1.50000 & 1.16089 \\ 3 & 1.15713 & 1.10340 \\ 4 & 1.09582 & 1.07582 \\ 5 & 1.06911 & 1.05969 \\ 6 & 1.05408 & 1.04913 \\ 7 & 1.04443 & 1.04170 \\ 8 & 1.03771 & 1.03619 \\ 9 & 1.03276 & 1.03195 \\ 10 & 1.02896 & 1.02858 \\ 20 & 1.01341 & 1.01381 \\ \end{array} \right)$$

Edit

More than inspired by @Svyatoslav's approach, let $$x_n=\sqrt{4 \sqrt{n}-2}+\sum_{k=1}^{12} \frac {a_k}{n^{(2k-1)/4}}$$

This choice of powers was done since the previous work showed that the coefficient of powers $n^{k/2}$ are all zero.

Using series, we have all coefficients functions of $a_1$ which is the solution of a polynomial of degree $13$ which shows only a real root very close to $\frac 12$ (more exactly $0.499973$).

As a result $$x_n \sim \sqrt{4 \sqrt{n}-2}+\frac 1{2\,n^{1/4}}-\frac 5{16\,n^{3/4}}+\frac 1{64\,n^{/4}}-\frac{301}{3072\,n^{7/4}}+$$ $$\frac{7}{4096\,n^{9/4}}-\frac{4265}{98304\,n^{11/4}}+\frac{33}{131072\,n^{13/4}}-\frac{6124957}{264241152\,n^{15/4}}+$$ $$\frac{715}{16777216\,n^{17/4}}-\frac{1099112215}{76101451776\,n^{19/4}}+\frac{4199}{536870912\,n^{21/4}}-$$ $$\frac{75614254421}{7653631721472\,n^{23/4}}+O\left(\frac{1}{n^{25/4}}\right)$$

Then, a few numbers $$\left( \begin{array}{ccc} n & \text{approximation} & \text{exact} & \text{|difference|} \\ 2 & 2.11586 & 2.12132 & 0.00546 \\ 3 & 2.44998 & 2.45465 & 0.00468 \\ 4 & 2.68563 & 2.68986 & 0.00423 \\ 5 & 2.87180 & 2.87574 & 0.00394 \\ 6 & 3.02753 & 3.03126 & 0.00373 \\ 7 & 3.16237 & 3.16594 & 0.00356 \\ 8 & 3.28189 & 3.28532 & 0.00343 \\ 9 & 3.38962 & 3.39294 & 0.00331 \\ 10 & 3.48796 & 3.49118 & 0.00322 \\ 20 & 4.18928 & 4.19195 & 0.00267 \\ \end{array} \right)$$