Solving recurrence sequence using an ODE

Question

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 ?

Answer 1

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:

Answer 2

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]:

Answer 3

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}, 
      g[q] /. Thread[
        C /@ Range[q - 1] -> data2[[q + 1 ;; 2 q - 1, q - 1]]]]

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}], 
  PlotStyle -> Red], plot, 
 PlotRange -> All]

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