Asymptotic analysis of $x_{n+1} = \frac{x_n}{n^2} + \frac{n^2}{x_n} + 2$

Question

Problem: Let $x_1 = 1$ and $x_{n+1} = \frac{x_n}{n^2} + \frac{n^2}{x_n} + 2, \ n\ge 1$. Find the third term in the asymptotic expansion of $x_n$.

I have posted it in MSE six months ago without solution for the third term https://math.stackexchange.com/questions/3801405/the-limit-and-asymptotic-analysis-of-a-n2-n-from-a-n1-fraca-nn.

We have $\lim_{n\to \infty} (x_n - n) = \frac{1}{2}$ (see [1]; I also give a solution with the help of computer in the link above). So the first two terms in the asymptotic expansion of $x_n$ are $x_n \sim n + \frac{1}{2}$.
Edit: In [1], the authors proved that $\frac{1}{4n-2} \le x_n - n - \frac{1}{2} \le \frac{2}{2n-3}$ for all $n\ge 3$.

For the third term, @Diger in MSE said $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$ (see @Diger's answer in the link above). However, I did some numerical experiment which does not support this result. I am not convinced of the numerical evidence due to finite precision arithmetic. I hope to prove or disprove it analytically.

Numerical Experiment: If $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$, then it should hold $16n(x_{2n} - 2n - \frac{1}{2}) \approx 5$ and $8(2n+1)(x_{2n+1} - (2n+1) - \frac{1}{2}) \approx 5$ for large $n$. When $n=1500$, I use Maple to get $16n(x_{2n} - 2n - \frac{1}{2}) \approx 4.368$ and $8(2n+1)(x_{2n+1} - (2n+1) - \frac{1}{2}) \approx 5.642$. When $n$ is larger (e.g., $n=10000$), the numerical result seems unreliable.

I ${\color{blue}{\textbf{GUESS}}}$ that $$x_{2n} \sim 2n + \frac{1}{2} + \frac{q_1}{2n},$$ $$x_{2n+1} \sim (2n+1) + \frac{1}{2} + \frac{q_2}{2n+1}$$ where $q_1 + q_2 = \frac{5}{4}$ and $q_1 \ne q_2$ (if $q_1 = q_2$, then it is $x_n \sim n + \frac{1}{2} + \frac{5}{8n}$). (Some numerical experiment shows $q_1 \approx \frac{61}{112}, q_2 \approx \frac{79}{112}$. But I am not convinced of it.)
Edit: I give more analysis for my guess as an answer.

Any comments and solutions are welcome and appreciated.

Reference

[1] Yuming Chen, Olaf Krafft and Martin Schaefer, “Variation of a Ukrainian Olympiad Problem: 10982”, The American Mathematical Monthly, Vol. 111, No. 7 (Aug. - Sep., 2004), pp. 631-632

Answer 1

Consider the substitutions \begin{equation*} x_n=n+1/2+y_n/n,\quad y_n=u_n+5/8. \end{equation*} Then $u_1=-9/8$ and \begin{equation*} u_{n+1}=f_n(u_n) \end{equation*} for $n\ge1$, where \begin{equation*} f_n(u):=\frac{-64 n^4 u-8 n^3 (4 u-13)+n^2 (56 u+115)+n (96 u+76)+4 (8 u+5)}{8 n^2 \left(8 n^2+4 n+8 u+5\right)}. \end{equation*} Define $c_n(u)$ by the identity \begin{equation*} f_n(u)=-u+\frac{13}{8n}+\frac{c_n(u)}{n^2}, \end{equation*} so that \begin{equation*} c_n(u)=\frac{n^2 \left(64 u^2+96 u+63\right)+n (11-8 u)+4 (8 u+5)}{8 \left(8 n^2+4 n+8 u+5\right)}. \end{equation*} Then for $n\ge1$ \begin{equation*} u_{n+1}+u_n=\frac{13}{8n}+\frac{c_n(u)}{n^2} \tag{1} \end{equation*} and for $n\ge2$ \begin{equation*} u_{n+1}=f_n(f_{n-1}(u_{n-1}))=u_{n-1}-\frac{13}{8n(n-1)}+\frac{c_n(u_n)}{n^2} -\frac{c_{n-1}(u_{n-1})}{(n-1)^2}. \tag{2} \end{equation*}

Note that \begin{equation*} u_{101}=-0.54\ldots,\quad u_{102}=0.56\ldots, \tag{3} \end{equation*} and \begin{equation*} 0\le c_n(u)\le3 \end{equation*} if $n\ge10$ and $u\in[-6/10,8/10]$. Therefore and because for natural $m\ge102$ we have \begin{equation*} \sum_{n=m}^\infty\Big(\frac{13}{8n(n-1)}+\frac3{(n-1)^2}\Big)<\frac5{m-2}\le0.05, \end{equation*} it follows from (2) and (3) by induction that for all $n\ge101$ we have $u_n\in[-6/10,8/10]$ and hence $0\le c_n(u_n)\le3$. So, again by (2), the sequences $(u_{2m})$ and $(u_{2m+1})$ are Cauchy-convergent and hence convergent. Moreover, by (1), $u_{n+1}+u_n\to0$.

Thus, indeed \begin{equation*} y_{n+1}+y_n\to5/4, \end{equation*} and the sequences $(y_{2m})$ and $(y_{2m+1})$ are convergent. (The limits of these two sequences can in principle be found numerically with any degree of accuracy -- controlled by (2), say.)

Answer 2

Something about my ${\color{blue}{\textbf{GUESS}}}$ for $q_1 = \frac{61}{112}, q_2 = \frac{79}{112}$:

Using the approach of @Iosif Pinelis, we consider the substitutions: for $n = 1, 2, \cdots$ \begin{align} x_{2n} &= 2n + \frac{1}{2} + \frac{\frac{61}{112} + u_n}{2n},\\ x_{2n+1} &= (2n+1) + \frac{1}{2} + \frac{\frac{79}{112} + v_n}{2n+1}. \end{align} We have $u_n + v_n \to 0$ as $n \to \infty$.

According to $x_{2n+1} = \frac{x_{2n}}{(2n)^2} + \frac{(2n)^2}{x_{2n}} + 2$ and $x_{2n+2} = \frac{x_{2n+1}}{(2n+1)^2} + \frac{(2n+1)^2}{x_{2n+1}} + 2$, we have \begin{align} v_n &= f(n, u_n), \\ u_{n+1} &= g(n, v_n) \end{align} where \begin{align} f(n, u) &= \frac{f_1(n,u)}{401408n^5 + 100352n^4 + (100352u+54656)n^3}, \\ f_1(n,u) &= -401408 n^5 u+(-100352 u+326144) n^4+(180096 u+223528) n^3 \\ &\quad +(150528 u+94528) n^2+(25088 u^2+52416 u+21106) n \\ &\quad +12544 u^2+13664 u+3721,\\ g(n, v) &= \frac{g_1(n,v)}{g_2(n,v)},\\ g_1(n, v) &= -401408 n^5 v+(-1103872 v+326144) n^4+(-1008000 v+916136) n^3 \\ &\quad +(-207424 v+1003708) n^2+(25088 v^2+182560 v+509856) n \\ &\quad +25088 v^2+91280 v+104375, \\ g_2(n, v) &= 401408 n^5+1103872 n^4+(100352 v+1275008) n^3\\ &\quad +(150528 v+758464) n^2 +(75264 v+228704) n+12544 v+27664. \end{align} We have $u_1 = 275/112$ and $$u_{n+1} = g(n, f(n, u_n)), \ n\ge 1.$$ I use Maple to get $u_{101} \approx 0.00342$, $u_{1001} \approx 0.000676$, $u_{10001} \approx 0.0004$.

Is it $\lim_{n\to \infty} u_n = 0$?

Answer 3

This is only a suggestion about how to get more terms of the asymptotic expansion.

Numerical experiments suggest that the sequences $x_{2n}$ and $x_{2n+1}$ have different asymptotic expression of type \begin{align*} x_{2n}&\sim 2n+\frac12+ \frac{p}{2n}\\ x_{2n+1}&\sim 2n+1+\frac12+\frac{q}{2n+1}. \end{align*}

We may assume that there are two power series $$f(x)=2+\frac{x}{2}+\sum_{k=2}^\infty a_k x^k, \quad g(x)=2+\frac{3x}{2}+\sum_{k=2}^\infty b_k x^k,$$ such that $$f(\tfrac{1}{n})=\frac{x_{2n}}{n},\qquad \frac{x_{2n+1}}{n}=g(\tfrac1n). $$ Then the recurrence gives us $$\tfrac14\tfrac{1}{n} f(\tfrac1n)+\frac{4}{\frac1nf(\tfrac1n)}+2=\frac{x_{2n}}{4n^2}+\frac{4n^2}{x_{2n}}+2=x_{2n+1}=\frac{g(\frac1n)}{\frac1n}$$ $$(n+1)f(\tfrac{1}{n+1})=x_{2n+2}=\frac{x_{2n+1}}{(2n+1)^2}+\frac{(2n+1)^2}{x_{2n+1}}+2=\frac{n g(\frac1n)}{(2n+1)^2}+\frac{(2n+1)^2}{n g(\frac1n)}+2.$$ If we put $z=\frac1n$ these means that the power series has to satisfy $$\frac{g(z)}{z}=\frac{z}{4}f(z)+\frac{4}{z f(z)}+2,\qquad \frac{z+1}{z} f(\tfrac{z}{z+1})=\frac{z g(z)}{(2+z)^2}+\frac{(2+z)^2}{zg(z)}+2.$$ If convergent power series satisfying these conditions exist, then $$x_{2n}=n f(\tfrac1n),\qquad x_{2n+1}=n g(n), $$ satisfies the recurrence and we will have finished.

Mathematica gives us $$f(x)=2+\frac{x}{2}+ a x^2+\frac{20a-3}{16}x^3+\frac{16a^2+74a-15}{64}x^4+\frac{1056a^2+1348 a-335}{1536}x^5+\cdots$$ $$g(x)=2+\frac{3x}{2}+\frac{-8a+5}{8}x^2-\frac{3(8a-3)}{32}x^3+\frac{32a^2-60a+19}{128}x^4+\frac{240a^2-266a+67}{768} x^5+\cdots$$

Different initial values probably leads to different values of $a$

With Mathematica we find $$a=0.27250624631159814997459302237509679293585\dots$$ by equating the first five terms of the expansion with $x_{10^7}/5\times10^6$. I have not make any other computations to try to get its approximation.

This number is not recognized by Mathematica as a simple algebraic number.