# Solving recursive equations with RSolve, with constraint and subscripts

I'm trying to solve this set of equations in generality, where $$n$$ can vary.

\begin{align*} P_{1}&=1-\sum_{i=2}^{n}P_{i}\\ P_{i}&=\frac{c_{i-1}P_{i-1}}{c_{i}+\mu_{i}}, \quad i=2,\dots,n-1\\ P_{n}&=\frac{c_{n-1}P_{n-1}}{\mu_{n}} \end{align*}

I've had luck solving specific special cases, for example when $$n=3$$, I've used:

Solve[
P1 == 1 - (P2 + P3) &&
P2 == (c1*P1)/(c2 + mu2) &&
P3 == (c2*P2)/(mu3),
{P1, P2, P3}
] // FullSimplify

{{P1 -> ((c2 + mu2) mu3)/(c1 c2 + (c1 + c2 + mu2) mu3),
P2 -> (c1 mu3)/(c1 c2 + (c1 + c2 + mu2) mu3),
P3 -> (c1 c2)/(c1 c2 + (c1 + c2 + mu2) mu3)}}


However, I'm struggling to use RSolve to do it in generality, I'm not sure if the issue has to do with my use of Indexed on the coefficients c and mu, or something else. My attempt is in the code block below. Many thanks in advance to any issues that can be pointed out!

eqns = {
P[nmax] == (Indexed[c, nmax - 1]*P[nmax - 1])/Indexed[mu, nmax],
P[n] == (Indexed[c, n - 1]*P[n - 1])/(Indexed[c, n] + Indexed[mu, n]),
P[1] == 1 - Sum[P[i], {i, 2, nmax}]
};

RSolve[eqns, P[n], n]


This is probably not the solution you are expecting, but note that this can be recast as a linear algebra problem; specifically, you are asking for the first column of a certain matrix inverse:

With[{n = 3},
LinearSolve[SparseArray[{Band[{1, 1}] -> 1, {1, k_} :> 1,
Band[{2, 1}] ->
Append[Table[-c[i - 1]/(c[i] + μ[i]), {i, 2, n - 1}],
-c[n - 1]/μ[n]]}, {n, n}], UnitVector[n, 1]]]


On a mathematical note, peering at the structure of the underlying matrix shows that it is the sum of a lower bidiagonal matrix and a rank-$$1$$ matrix, which allows you to use the Sherman-Morrison-Woodbury formula.

Letting $$\mathbf A=\mathbf B+\mathbf e_1(\mathbf e-\mathbf e_1)^\top$$, where $$\mathbf B$$ is the lower bidiagonal part, $$\mathbf e$$ is the vector with components all equal to $$1$$, and $$\mathbf e_1$$ is the first unit vector, you have the relation

\begin{align*} \mathbf A^{-1}\mathbf e_1&=\mathbf B^{-1}\mathbf e_1-\frac{(\mathbf B^{-1}\mathbf e_1)(\mathbf e-\mathbf e_1)^\top(\mathbf B^{-1}\mathbf e_1)}{1+(\mathbf e-\mathbf e_1)^\top\mathbf B^{-1}\mathbf e_1}\\ &=\frac1{1+(\mathbf e-\mathbf e_1)^\top\mathbf B^{-1}\mathbf e_1}\mathbf B^{-1}\mathbf e_1\\ &=\frac{\mathbf B^{-1}\mathbf e_1}{\mathbf e^\top\mathbf B^{-1}\mathbf e_1} \end{align*}

Note that $$\mathbf B$$ is a lower triangular matrix with a special structure (specifically, a semiseparable matrix, see e.g. the book of Vandebril et al.); you might be able to perform a further simplification afterwards.

In Mathematica, this means that the following code should give the same result:

With[{n = 3},
Normalize[LinearSolve[SparseArray[{Band[{1, 1}] -> 1,
Band[{2, 1}] ->
Append[Table[-c[i - 1]/(c[i] + μ[i]), {i, 2, n - 1}],
-c[n - 1]/μ[n]]}, {n, n}],
UnitVector[n, 1]], Total]]


Using the special structure of the bidiagonal inverse, as noted earlier, this means that the following code should also be equivalent:

With[{n = 3},
Normalize[FoldList[Times, 1, Append[Table[c[i - 1]/(c[i] + μ[i]), {i, 2, n - 1}],
c[n - 1]/μ[n]]], Total]]

• thank you, this is fascinating and much more elegant when considered in this way. – slwu89 Mar 7 '19 at 6:16