Convergence or divergence of $a_{n+1}=a_n +\frac{a_{n-1}}{(n+1)^2}$

Question

If $a_1=a_2=1$ and $a_{n+1}=a_n +\frac{a_{n-1}}{(n+1)^2}$ How to prove convergence of the sequence and its limit or divergence?

It is easy to see that the sequence is always positive and by that one can conclude it is strictly increasing sequence but proving an upper bound doesn't seem to be easy, so maybe it diverge?

I wrote a simple python code and found that all $a_n =1.4552363391830$ where $n = 200,000,000$ this might suggest that the sequence is bounded but it is worth noting that this sequence might have a very slow convergence or a very slow divergence because $a_n=1.4552\color{red}{291356494609}$ where $n=200,000$ so the sequence is able to gain about $10^{-5}$ from $2\times 10^5$ to $2\times 10^8$.

import math
x=1
y=1
for i in range(3, 1_000_000_000):
    x=y+x/((2*i-3)**2)
    y=x+y/((2*i-2)**2)    
print(x)
print(y)

$$a_{2*10^9} =1.4552364010671075$$

If this sequence diverges then how many terms one need for $a_n >x$ for any $x\in \mathbb{R}^+$? If this sequence converges can a closed form for a the limit exist ?

Answer 1

$$a_{n+1} = a_n + \frac{a_{n-1}}{(n+1)^2} < a_n\left(1 + \frac{1}{(n+1)^2}\right)$$ $$\frac{a_{n+1}}{a_{n}} < 1 + \frac{1}{(n+1)^2}$$ $$\frac{a_{n}}{a_{1}} < \prod_{k=2}^n\left (1 + \frac{1}{k^2}\right) = e^{\sum_{k=2}^n \ln(1 + \frac{1}{k^2})} < e^{\sum_{k=2}^n \frac{1}{k^2}} \to e^{\frac{\pi^2}{6}-1}$$ $$\lim_{n \to +\infty} a_{n} < e^{\frac{\pi^2}{6}-1} \approx 1.905861$$

Answer 2

Let $S$ be a real number greater than one (to be determined later) and assume that $a_n \le S$ for $1 \le n \le N$. Then $$ a_{N+1} = 1 + \sum_{n=2}^N (a_{n+1}-a_n) = 1 + \sum_{n=2}^N \frac{a_{n-1}}{(n+1)^2} \\ \le 1 + \sum_{n=2}^\infty \frac{S}{(n+1)^2} = 1 + S \left( \frac{\pi^2}{6} - 1 - \frac 14\right) \, . $$ Now we choose $S$ such that the right-hand side of the last inequality is equal to $S$, that is $$ S = \frac{1}{9/4 - \pi^2/6} \approx 1.652712 \, . $$ Then the above calculation shows that $a_n \le S$ for $1 \le n \le N$ implies $a_{N+1} \le S$, so that $a_n \le S$ for all $N$ by induction.

This proves that $(a_n)$ is convergent, and $$ \lim_{n\to \infty} a_n \le \frac{1}{9/4 - \pi^2/6} \approx 1.652712 \, . $$

---

The bound can be improved if we compute some initial values. For example with $$ a_1 = a_2 =1 , a_3 = \frac{10}{9}, a_4 = \frac{169}{144} $$ we get that $a_n \le S$ for $1 \le n \le N$ implies $$ a_{N+1} \le \frac{169}{144} + S \left( \frac{\pi^2}{6} - 1 - \frac 14 - \frac 19 - \frac{1}{16}\right) $$ and that gives the upper bound $$ S = \frac{169}{349-24 \pi^2} \approx 1.507186 \ . $$

Answer 3

This is a note on high-precision computation of the limit $$a_\infty=\lim\limits_{n\to\infty}a_n\approx1.45523641194065421249837906730024813\cdots$$

---

First, let $a_n=a_\infty+b_n$. Then $n^2(b_n-b_{n-1})=a_\infty+b_{n-2}$, hence $$ a_\infty=\lim_{n\to\infty}n^2(b_n-b_{n-1})=-\lim_{n\to\infty}\frac{b_n-b_{n-1}}{\frac1n-\frac1{n-1}}, $$ and the Stolz–Cesàro theorem yields $\lim\limits_{n\to\infty}nb_n=-a_\infty$. This alone already gives a considerable speedup: $a_n^{(1)}=\frac{na_n}{n-1}$ converges to $a_\infty$ much faster than $a_n$ itself (see the table below).

---

But this process can be continued. The result is the asymptotic expansion $$ a_n\asymp a_\infty\sum_{k=0}^{(\infty)}c_k n^{-k}\qquad(n\to\infty) $$ with the coefficients $c_k$ given by $c_0=1$, $c_1=-1$, $c_2=1$ and, for $k>1$, $$ c_{k+1}=-\frac1{k+1}\sum_{j=1}^k\left[\binom{k+1}{j-1}+2^{k-j}\binom{k-1}{j-1}\right]c_j. $$ So, the sequence $(c_k)$ begins with $$ 1,-1,1,-\frac13,-\frac16,-\frac16,-\frac19,\frac{19}{63},\frac{227}{168},\frac{13499}{4536},\frac{16183}{4536},-\frac{220859}{49896},-\frac{14209739}{299376},\dots\newcommand{\ntail}[1]{\color{\LightBlue}{#1\cdots}} $$

---

Now we can compute $a_n^{(m)}=a_n/\sum_{k=0}^m c_k n^{-k}$ instead of $a_n$:

$n$	$a_n$	$a_n^{(1)}$	$a_n^{(2)}$	$a_n^{(20)}$
$10^2$	$1.4\ntail{408290}$	$1.455\ntail{38291306447}$	$1.45523\ntail{5919537251482034}$	$38$ digits
$10^3$	$1.45\ntail{37826}$	$1.45523\ntail{786814795}$	$1.455236411\ntail{454847305315}$	$58$ digits
$10^4$	$1.455\ntail{0909}$	$1.4552364\ntail{2649398}$	$1.455236411940\ntail{169060927}$	$79$ digits
$10^5$	$1.4552\ntail{218}$	$1.45523641\ntail{208617}$	$1.45523641194065\ntail{3727412}$	$99$ digits
$10^6$	$1.45523\ntail{49}$	$1.45523641194\ntail{210}$	$1.455236411940654212\ntail{013}$	$121$ digits
$10^7$	$1.455236\ntail{2}$	$1.4552364119406\ntail{6}$	$1.45523641194065421249\ntail{7}$	$142$ digits

---

Another approach uses generating functions.

Let $g(z):=\sum_{n=1}^\infty a_n z^n$ for $|z|<1$; then we have $a_\infty=\lim\limits_{z\to1^-}(1-z)g(z)$.

Let $w(z):=(1-z)g(z)=\sum_{n=1}^\infty c_n z^n$ where $c_1=1$ and $c_n=a_n-a_{n-1}$ for $n>1$. The given recurrence for $a_n$ is now rewritten as $(n+1)^2 c_{n+1}=\sum_{k=1}^{n-1}c_k$, which holds for $n>0$ (under the agreement that the sum of zero terms is zero). Multiply by $z^n$ and sum over all $n>0$; since $$ \sum_{n=1}^\infty(n+1)^2 c_{n+1}z^n=\frac{d}{dz}\left(z\frac{d}{dz}\sum_{n=1}^\infty c_{n+1}z^{n+1}\right)=\big(zw'(z)\big)'-1, $$ $$ \sum_{n=1}^\infty z^n\sum_{k=1}^{n-1}c_k=\sum_{k=1}^\infty c_k\sum_{n=k+1}^\infty z^n=\sum_{k=1}^\infty\frac{c_k z^{k+1}}{1-z}=\frac{zw(z)}{1-z}, $$ we obtain $\color{blue}{z(1-z)w''(z)+(1-z)w'(z)-zw(z)=1-z}$, and seek for $a_\infty=w(1^-)$.

This can be solved using modern numerical methods for ODEs (sinc-collocation with double-exponential transformations), which has a chance to outperform the preceding approach.

Answer 4

The sequence $a_n$ is strictly increasing. So, $a_{n+1}>a_n,~\forall n>1$.

Now, $$a_{n+1}=a_n+\frac{a_{n-1}}{(n+1)^2}<a_n+\frac{a_n}{(n+1)^2}=a_n\left(1+\frac{1}{(n+1)^2}\right)\tag{1}$$

Taking natural logarithms on both sides and applying $\ln(1+x)<x,\forall x>0$, we get $$\ln a_{n+1}<\ln a_n+\ln\left(1+\frac{1}{(n+1)^2}\right)<\ln a_n+\frac{1}{(n+1)^2}\tag{2}$$

Taking summation on both sides, $$\sum_{n=1}^N\left(\ln a_{n+1}-\ln a_n\right)<\sum_{n=1}^N\frac{1}{(n+1)^2}<\sum_{n=1}^\infty\frac{1}{(n+1)^2}=\zeta(2)-1=\frac{\pi^2}6-1$$ The left-hand side telescopes and the upper bound of $(\ln a_{N+1}-\ln a_1)$ or $\ln a_{N+1}$ is $\frac{\pi^2}6-1$. That is, an upper bound of $a_{N+1}$ is $\exp\left(\frac{\pi^2}6-1\right)\approx1.90586$. Hence, the sequence defined by the given recurrence relation converges and $a_n\to\ell<e^{\zeta(2)-1}$ as $n\to\infty$, which is consistent with the result by your Python code. $\blacksquare$