Finding Steady States that Meet Additional Constraints

Question

I've written a function that returns the steady state given a regular, column-stochastic matrix. I want to use it to solve a larger problem, finding steady states that satisfy certain conditions.

The function seems to work fine on its own. But not as part of the larger problem.

Here's the function definition:

SteadyState[M_] := 
 First[Eigenvectors[M]] / Total[First[Eigenvectors[M]]]

Used on its own it seems to work well, e.g.:

SteadyState[
  {{1/2, 1/2},
   {1/2, 1/2}}
]

returns {1/2, 1/2} as it should.

But when I try to use it inside of Solve, it doesn't behave as expected. For example:

Solve[
 SteadyState[
    {{c11, c12},
     {c21, c22}}
    ] == {q1, q2} &&
  q1 >= 0 && q2 >= 0 &&
  q1 + q2 == 1 &&
  c11 >= 0 && c12 >= 0 && c21 >= 0 && c22 >= 0 &&
  c11 + c21 == 1 && 
  c12 + c22 == 1 &&
  c11 == 1/2 && c12 == 1/2 && q1 == 1/2 && q2 == 1/2,
 {q1, q2, c11, c12, c21, c22},
 Reals
 ]

Returns the empty set of solutions {}, even though I've effectively handed it the solution in the last two constraints!

Could this be because my SteadyState function doesn't check that the input is regular and column stochastic? If so, why? And what would be a good way to rectify that?

Answer 1

I should have seen right away, for this purpose it's sufficient—and way more sensible—to just put the steady-state condition into Solve by hand, as one of the constraints.

To simplify the problem's expression I first defined some helpers:

ColumnStochastic[M_] := AllTrue[
   Table[i, {i, Length[M]}],
   Total[M[[All, #]]] == 1 &
   ] && AllTrue[Flatten[M], # >= 0 &]

pvec[v_] := AllTrue[v, # >= 0 &] && Total[v] == 1

The following code then runs no problem:

A = {{a11, a12}, {a21, a22}};
q = {q1, q2};

FindInstance[
  A.q == q &&
    ColumnStochastic[A] && pvec[q],
  Join[q, Flatten[A]]
]

It finds the solution:

{{q1 -> 3/4, q2 -> 1/4, a11 -> 27/32, a12 -> 15/32, a21 -> 5/32, a22 -> 17/32}}

Using Solve or Reduce instead of FindInstance also works now to get the general class of solutions.

To find steady-states that meet additional constraints (as alluded in the OP), I can now just && them in alongside the others.