Solving a system of Linear Equations with constraints. Using Symbolic math.

I have the following system of equations with a constraint:

n equations with n unknowns

Here are all known. (Its a known n*n matrix of values).

The following is my try at a solution. How can I correct it to get the 's

p_ij = [0.5, 0.5; 0.2, 0.8];

length_pi_i = size(p_ij);
length_pi_i = length_pi_i(1);

syms pi [1 length(p_ij)]
assume(sum(pi) == 1);


zero_vector = zeros(length(pi), 1);
eqn = (pi' - (pi * p_ij)' == zero_vector);


M = solve(eqn, pi);
disp(M);

Answer 1

[Tanmay Das]

Hi,

I tried reproducing your steps in MATLAB R2021b release and received the following output in the command window:

pi1: 2/7

pi2: 5/7

I tried verifying the answer manually and the answer seems correct to me.

In case you want to learn more about symbolic functions, you may look into 'syms','assume', and 'solve' function documentation.

Answer 2

[John D'Errico]

This is just an eigenvalue problem. (Yes, I know that is probably beyond the scope of the question that was asked. But this is just a basic Markov chain problem. And we should recognize the matrix P_ij as a 2x2 Markov transition matrix.)

p_ij = [0.5, 0.5; 0.2, 0.8]

p_ij = 2×2

0.5000 0.5000 0.2000 0.8000

But first, we can ask if a solution of the form requested exists for this problem? That is, does the problem

X*P_ij = X

Where X is a vector of length 2, have a solution? For that to be the case, then we must have 1 as an eigenvalue of this matrix.

format long g
[V,D] = eig(p_ij.')

V = 2×2

-0.707106781186547 -0.371390676354104 0.707106781186547 -0.928476690885259

D = 2×2

0.3 0 0 1

So in fact, it is true that 1 is an eigenvalue. We see that because we see that D(2,2) == 1. The corresponding vector, thus V(:,2) is the solution.

pi_vec = V(:,2)'

pi_vec = 1×2

-0.371390676354104 -0.928476690885259

But we should also know that an eigenvector can be scaled arbitrarily. We can multiply it by any constant, and it is still an equally valid eigenvector. So we can now just rescale pi_vec to sum to 1 as needed. And by switching to format rat, we find the desired long term probabilities for each of the two states of this Markov process.

format rat
pi_vec = pi_vec/sum(pi_vec)

pi_vec = 
       2/7            5/7     

There was not even a need to use solve. Just eig was sufficient here. Could it have been done using syms? Of course, though I refuse to name a variable pi, as that tramples on the terribly useful number pi.

syms X [1,2]
T = sym(p_ij) %% the transition matrix, in symbolic form.

T = 

solve(X*T == X,sum(X) == 1,X)

ans = struct with fields:
    X1: 2/7
    X2: 5/7

Again, we see the expected solution drops out.