# Solving recurrence sequence using an ODE

I'm trying to solve the recurrence

a[n] == 0 for n < l
a[l] == 1
a[n] == (l + 2(n - l) a[n - l] + (n - l)(n - 1) a[n - 1])/n(n - l + 1) for n >= l


using the ODE

l = 2;
DSolve[(x^2 y''[x] + (2 - l) x y[x]) (x - 1) + 2 x^(l + 1) y'[x] + (l x^l)/(1 - x) == 0,
y[x], x]


For l = 2 it works, but not for l > 2. Is it possible to do it? Maybe there is another way to solve this recurrence ?

-
@m_goldberg The previous version was useful for experimenting with RSolve[] – Dr. belisarius Nov 11 '13 at 16:42
Which is the real problem you need to solve: the recurrence or the ODE? – bill s Nov 14 '13 at 13:31
@bills I need to solve the recurrence. – Филипп Цветков Nov 14 '13 at 13:50
You mean you want to get a analytic formula of a[n] for any l? – xzczd Nov 16 '13 at 3:12
Could you please post the solution for a[n] for l=2 ? I don't know how to link the second equation to the first. – b.gatessucks Nov 17 '13 at 11:16

After a change of Variable

s[k_]:= s[k]= RSolve[{a[p + m] == ((1 + p) (m + 2 p a[p] + p (-1 + m + p) a[-1 + m + p]))/(m + p) /.
m -> k,  Sequence @@ Array[a@# == 0 &, k - 1, 0]}, a, p]


Which give DifferenceRoot solutions. So:

Grid[Table[a[n] /. s[m][[1]], {n, 1, 5}, {m, 1, 5}]] /. C[1] -> 0


Gives the same table as in @ubpdqn's answer, except that for l==1 you get the solution expressed as some combinations of Bessel functions:

-
@belisarius... compact...reassuring for me. – ubpdqn Nov 19 '13 at 8:13
I think that's the best one can do (+1) – gpap Nov 19 '13 at 13:02

This is not an analytic solution of the recurrence. I generate a table based on the recurrence relation (if I interpret it correctly). I have change l to m as difficult to discriminate from one. If the values are wrong I have obviously made an error but perhaps can motivate other approaches.

f[m_, n_, mat_] :=
1/n(1 - m + n) (m + (-1 + n) (-m + n) mat[[m, -1 + n]] +
2 (-m + n) mat[[m, -m + n]])
g[n_] := Module[
{sa, pos},
sa = Normal@
SparseArray[{{i_, i_} -> 1, {i_, j_} /; i > j -> 0}, {n, n}, "x"];
pos = Position[sa, "x"];
Fold[ReplacePart[#1, #2 -> f[#2[[1]], #2[[2]], #1]] &, sa, pos]

]


Tabulating g[5]:

-

Are you not satisfied with difference root objects? I think that @belisarius has given the most general solution. If you want to give your sequence values and see how it behaves there's no reason to use RSolve. You can define the n-th term recursively and build the sequence from the bottom up (in what follows, I changed l to q):

ClearAll[a, g, data, data2, plot, an6];
a[n_] := If[n >= q,
((1 + n - q) (q + (-1 + n) (n - q) a[-1 + n] + 2 (n - q) a[n - q]))/
n,
0]


Now you can generate a few terms and watch a grow:

data = Table[Table[{q, k, Log[a[k]]}, {q, 1, k}], {k, 1, 16}];


where I take the logarithm because it grows too fast and I want to plot it. Here's how (the logarithm of) your sequence looks for the first few n, q:

plot = ListPointPlot3D[Flatten[data, 1]]


If you are after solutions for large n and small qs the following might be a bit quicker than @belisaruis's:

    g[q_] := b[n] /.
RSolve[{b[n] == ((1 + n - q) (q + (-1 + n) (n - q) b[-1 + n] + 2 (n - q) b[n - q]))/n,
b[q] == 1},
b[n], n]


but for every g[n] it needs n-1 additional values to determine an equal number of constants. So

g[1]


evaluates to

(*out*){(1/(n! (2 + n)!))
2 (BesselI[2, 2] n! (2 + n)! -
HypergeometricPFQ[{1}, {1 + n, 3 + n}, 1]) Pochhammer[
1, -1 + n] Pochhammer[3, -1 + n]}


but

g[2]


evaluates to a difference root object with one constant. The point here is that you have another constant to determine for every increment of q. So I can use my a[n] above to determine the constants at every step:

data2 = Table[Table[a[k], {q, 1, k}], {k, 1, 17}];


which gives a table like the ones described in other answers. Now, say

an6 = First@With[{q = 6},

gives you a difference root object for q=6 and any n$\geq$q. The last bit is a set of rules that replaces the constants based on the table of values data2 above. This is exactly where the problem lies: the higher the q the more constants you need and my function a[n] gets quite slow for larger n, q so this essentially transposes the problem with speed to generating values to replace the constants. As I said, I don't think you can do anything more general than what @belisarius did. But if you care for the first few qs this might be more viable as you have a difference root object for any n given a q so in my example for q=6 above, you can go to any n:
Show[ListPointPlot3D[Table[{6, n, Log[an6]}, {n, 1, 35}],

Using data2 you can do the same thing for q between 1 and 9 and this will definitely be fast the more n grows.