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:

  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]]
  • $\begingroup$ thank you, this is fascinating and much more elegant when considered in this way. $\endgroup$ – slwu89 Mar 7 '19 at 6:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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