Solving recursive equations with RSolve, with constraint and subscripts

Question

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]

Answer 1

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