Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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 ?

share|improve this question
@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:

Mathematica graphics

share|improve this answer
@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]:

enter image description here

share|improve this answer

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]))/

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

Mathematica graphics

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


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



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]

Mathematica graphics

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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.