Find the formula for the sequence $(a_n)$ that satisfies the recurrence relation $a_n=(n+7)a_{n-1}+n^2$ with $a_0=1$

Question

Find the formula for the sequence $(a_n)$ that satisfies the recurrence relation $a_n=(n+7)a_{n-1}+n^2$ with the initial condition $a_0=1$.

This is a non-linear nonhomogeneous recurrence relation. What I can think that may help solve this problem is using the idea of generating function.
Let $G(x)=\sum_{n=0}^{\infty}a_nx^n$. Then $$\begin{aligned} G(x)&=a_0+\sum_{n=1}^{\infty}a_nx^n\\ &=1+\sum_{n=1}^{\infty}\left((n+7)a_{n-1}+n^2\right)x^n\\ &=1+\sum_{n=1}^{\infty}\left((n+7)a_{n-1}x^n\right)+\sum_{n=1}^{\infty}(n^2x^n) \end{aligned}$$ Then I was stuck here, because I couldn't change the two summations into forms of $G(x)$.
Anyone has brillant ideas?

Answer 1

Following @GTonyJacobs remark about the homogeneous sequence, let us define (this is a discrete method of variation of the parameter)

$${a}_{n} = \left(n+7\right) ! \ {b}_{n}$$

We have

$$\left(n+7\right) ! \ {b}_{n} = \left(n+7\right) \left(n+6\right) ! \ {b}_{n-1}+{n}^{2}$$

hence

$${b}_{n}-{b}_{n-1} = \frac{{n}^{2}}{\left(n+7\right) ! \ }$$

It follows that

$${b}_{n} = {b}_{0}+\sum _{k = 1}^{n} \left({b}_{k}-{b}_{k-1}\right) = \frac{1}{7 ! \ }+\sum _{k = 1}^{n} \frac{{k}^{2}}{\left(k+7\right) ! \ }$$

hence

$${a}_{n} = \left(n+7\right) ! \ \left(\frac{1}{7 ! \ }+\sum _{k = 1}^{n} \frac{{k}^{2}}{\left(k+7\right) ! \ }\right)$$

Answer 2

$x^8G(x) =\sum_{n=0}^{\infty}a_nx^{n+8}$ so $(x^8G(x))' =\sum_{n=0}^{\infty}(n+8)a_nx^{n+7} =\sum_{n=1}^{\infty}(n+7)a_{n-1}x^{n+6} =x^6\sum_{n=1}^{\infty}(n+7)a_{n-1}x^{n} $.

From this you can set up a differential equation for $G(x)$.

Answer 3

Since the solution of $$ b_{n} = (n+7) b_{n-1},\quad b_0=1 $$ is clearly given by $b_n=\frac{(n+7)!}{7!}$, it looks like a good idea to enforce the substitution $a_n=\frac{(n+7)!}{7!}A_n$, leading to $A_0=1$ and

$$ A_n - A_{n-1} = \frac{7!n^2}{(7+n)!}.\tag{A} $$ By summing both sides of $(A)$ over $n=1,2,\ldots,N$ we get $$ A_N = 1+\sum_{n=1}^{N}\frac{7!n^2}{(7+n)!}\tag{B} $$ and $$ a_N = \frac{(N+7)!}{7!}+\sum_{n=1}^{N}\frac{(N+7)!n^2}{(n+7)!}.\tag{C} $$ In particular, for large values of $N$ we have $$ a_N \approx (N+7)!\cdot\left(37e-\frac{21121}{210}\right). \tag{D}$$

Answer 4

Set $b_n=a_n+n-5$ to eliminate the term in $n^2$ and we get $b_n=(n+7)b_{n-1}+37$

Now set $b_n=(n+7)!c_n$ to get the difference formula $c_n-c_{n-1}=\dfrac{37}{(n+7)!}$

Sum the telescoping series and $c_n=c_0+37\sum\limits_{k=8}^{n+7}\dfrac 1{k!}$

With $a_0=1$ then $b_0=-4$ and $c_0=-\dfrac 4{7!}$ and after isolating the partial sum from $k=\{0..7\}$

$a_n=(n+7)!\left(37\sum\limits_{k=0}^{n+7}\dfrac 1{k!}-\dfrac{21121}{210}\right)+5-n$

We can eventually replace by $\sum\limits_{k=0}^{n+7}\dfrac 1{k!}=e\,\Gamma(n+8,1)$ the incomplete gamma function, explaining the equivalent given by Jack D'Aurizio.

Answer 5

This example is quite interesting because the generating function method leads to a highly singular integral. But let us analyze it.

\begin{eqnarray} G(x) &= & a_0 + \sum\limits_{n=1}^\infty a_n x^n \\ &=& a_0 + \sum\limits_{n=1}^\infty \left( (n+7) a_{n-1} + n^2 \right) x^n \\ &=& a_0 + \sum\limits_{n=0}^\infty (n+8) a_n x^{n+1} + \sum\limits_{n=0}^\infty n^2 x^n \\ &=& a_0 + x^2 d_x G(x) + 8 x G(x) + \frac{x(1+x)}{(1-x)^3} \end{eqnarray}

Therefore we have:

\begin{equation} d_x G(x) + \left(\frac{8}{x} - \frac{1}{x^2} \right) G(x) = - \frac{a_0}{x^2} - \frac{(1+x)}{x(1-x)^3} \end{equation}

The solution to the homogeneous ODE is straightforward and reads $x^{-8} e^{-1/x} $. Therefore the special solution to the in-homogeneous ODE reads:

\begin{equation} G(x) = x^{-8} e^{-1/x} \cdot \int x^8 e^{1/x} \left( - \frac{a_0}{x^2} - \frac{(1+x)}{x(1-x)^3}\right) dx \end{equation}

Now, the antiderivative is obtained by substituting for $1/x$ and then reducing the rational function in partial fractions and then integrating by parts. But this is tedious and therefore we used Mathematica to do those calculations. The result is the following:

\begin{equation} G(x) = \frac{(\text{a0}-506905) e^{-1/x} \text{Ei}\left(\frac{1}{x}\right)}{5040 x^8}+\frac{37 e^{1-\frac{1}{x}} \text{Ei}\left(\frac{1}{x}-1\right)}{x^8}+\frac{\frac{35 \left(24 x^7+72 x^6+138 x^5+298 x^4+825 x^3+3827 x^2-14483 x+9155\right)}{(x-1)^2}-\text{a0} \left(720 x^6+120 x^5+24 x^4+6 x^3+2 x^2+x+1\right)}{5040 x^7} \end{equation}

Now, we need to invert. This requires some knowledge about the behavior of the exponential integral function about plus infinity. We leave the full analysis of that for later. Here we only verify with Mathematica that the result is correct.

In[217]:= 
G[x_] = Collect[(Integrate[
     x^8 Exp[1/x] (-a0/x^2 - (1 + x)/(x (1 - x)^3)), 
     x] x^(-8) Exp[-1/x]), ExpIntegralEi[_], Simplify]
(D[G[x], x] + (8/x - 1/x^2) G[x]) // Simplify
Normal[Series[G[x], {x, 0, 3}, Assumptions -> x > 0]]

Out[217]= (-a0 (1 + x + 2 x^2 + 6 x^3 + 24 x^4 + 120 x^5 + 
     720 x^6) + (
  35 (9155 - 14483 x + 3827 x^2 + 825 x^3 + 298 x^4 + 138 x^5 + 
     72 x^6 + 24 x^7))/(-1 + x)^2)/(5040 x^7) + (
 37 E^(1 - 1/x)
   ExpIntegralEi[-1 + 1/x])/x^8 + ((-506905 + a0) E^(-1/x)
   ExpIntegralEi[1/x])/(5040 x^8)

Out[218]= (-a0 + (x (1 + x))/(-1 + x)^3)/x^2

Out[219]= a0 + (1 + 8 a0) x + (13 + 72 a0) x^2 + (139 + 720 a0) x^3