# Solve a set of recurrence relations with a constraint

I am trying to solve a set of recurrence equations as given below:

$p_{1, n}*(l + n*u1) = p_{1, n-1}*l + p_{1, n+1}*(n+1)*u1 + p_{2, n}*u2$,

$p_{2, n}*(l + n*u1 + u2) = p_{2, n-1}*l + p_{2, n+1}*(n+1)*u1,$

$p_{1, 1}*(l + n) = p_{2, 1}*u2 + p_{1, 2}*u1,$

$p_{2,C-1}*((C - 1)*u1 + u2) = p_{2, C-2}*l$,

$p_{1, C}*C*u1 = p_{1, C-1}*l$,

$\sum_{i=1}^C p_{1, i} + \sum_{i=1}^C p_{1, i} = 1$

I have used the following command :

RSolve[{p[1, n]*(l + n*u1) ==
p[1, n - 1]*l + p[1, n + 1]*(n + 1)*u1 + p[2, n]*u2,
p[2, n]*(l + n*u1 + u2) == p[2, n - 1]*l + p[2, n + 1]*(n + 1)*u1,
p[1, 1]*(l + n) == p[2, 1]*u2 + p[1, 2]*u1,
p[2, C - 1]*((C - 1)*u1 + u2) == p[2, C - 2]*l,
p[1, C]*C*u1 == p[1, C - 1]*l,
Sum[p[1, i], {i, 1, C}] + Sum[p[2, i], {i, 1, C - 1}] == 1}, {p[1,
n], p[2, n]}, n]


However, the output obtained is the same as the input, with the terms in the equation rearranged. What is wrong with what I am doing?

• It merely means RSolve[] doesn't know how to solve your system. Feb 20 '16 at 7:19
• @J.M. Thank you. What command can be used to solve such an equation? Feb 20 '16 at 7:20
• Are you even certain that there is a closed form solution to this? Most partial difference equations don't. Feb 20 '16 at 7:23
• @J.M. I have come across this system while solving for the stationary distribution of a 2D markov chain. There ought to be a solution. Feb 20 '16 at 7:24

Edit: Simplified Code

This problem certainly is solvable, but not with RSolve, which is not meant for such cases. (Also, n should not occur in the boundary conditions. I presume that it should be c and replace it accordingly.) Specify c and use Solve. For instance,

c = 4;
eqs = Join[
Table[p[1, n]*(l + n*u1) == p[1, n - 1]*l + p[1, n + 1]*(n + 1)*u1 + p[2, n]*u2,
{n, 2, c - 1}],
Table[p[2, n]*(l + n*u1 + u2) == p[2, n - 1]*l + p[2, n + 1]*(n + 1)*u1, {n, 2, c - 2}],
{p[1, 1] == norm, p[1, 1]*(l + c) == p[2, 1]*u2 + p[1, 2]*u1,
p[2, c - 1]*((c - 1)*u1 + u2) == p[2, c - 2]*l, p[1, c]*c*u1 == p[1, c - 1]*l}];
var = Join[Table[p[1, n], {n, c}], Table[p[2, n], {n, c - 1}]];
sol = Flatten@Solve[eqs, var];
gcd = PolynomialGCD @@ Denominator[Rest[var] /. sol];
sol = Simplify[sol /. norm -> norm*gcd]
rn = norm -> Simplify[norm/Total[var /. sol]]


The solution then is

{p[1, 1] -> norm (l^2 + 3 l (u1 + u2) + 2 (6 u1^2 + 5 u1 u2 + u2^2)),
p[1, 2] -> (l norm (l^2 + 6 u1^2 + u2 (4 + u2) + l (4 + 3 u1 + 2 u2) +
u1 (12 + 5 u2)))/u1,
p[1, 3] -> (l^2 norm (l^2 + 6 u1^2 + 5 u1 (4 + u2) + u2 (4 + u2) +
2 l (2 + 2 u1 + u2)))/(3 u1^2),
p[1, 4] -> (l^3 norm (l^2 + 6 u1^2 + 5 u1 (4 + u2) + u2 (4 + u2) +
2 l (2 + 2 u1 + u2)))/(12 u1^3),
p[2, 1] -> ((8 + l) norm (6 u1^2 + 5 u1 u2 + u2 (l + u2)))/u2,
p[2, 2] -> (l (8 + l) norm (3 u1 + u2))/u2,
p[2, 3] -> (l^2 (8 + l) norm)/u2}

norm -> (12 u1^3 u2)/(l^5 u2 + 2 l^4 u2 (2 + 4 u1 + u2) +
24 u1^3 (4 + u2) (6 u1^2 + 5 u1 u2 + u2^2) +
12 l u1^2 (6 u1^3 + u2^2 (4 + u2) + 2 u1^2 (12 + 7 u2) + u1 u2 (28 + 9 u2)) +
4 l^2 u1 (9 u1^3 + 24 u1^2 (1 + u2) + u2^2 (4 + u2) + u1 u2 (32 + 11 u2)) +
l^3 (12 u1^3 + 34 u1^2 u2 + u2^2 (4 + u2) + u1 u2 (36 + 13 u2)))


In symbolic form, this is cumbersome, even for c = 4. It simplifies enormously, if {u1, u2, l} are given numerical values.

• Thank you! Would it be possible to get a solution for general C? Feb 21 '16 at 15:38
• @PriyaP I can see patterns emerging as I increase c, but I do not know a formal way to derive a general solution. I am glad that you found what I was able to derive to be useful. Feb 22 '16 at 0:28