Solve linear equation in matrix form

I need to find pi^(0) from this linear equation

pi^(0) is a row vector, e is column vector of 1

where A_0, R, B, A, and C are square matrix.

Matrix A_0, B, A, and C are given.

R can be found by using successive substitutions method

After i found R, i use this code to find pi^(0) (q is the same as pi^(0))

q=sym('q',[1 mm]) %row vector

I=eye(mm);

e=ones(mm,1); %column vector of 1

a=q*(inv(I-R))*e %size 1 1

aa=a-1

j=solve(aa,sym('q1')) %j is the same as q1=

w=q*(A_0+R*B) %size 1 row mm column

k=subs(w,sym('q1'),j) %subs q1 with j

[C,S]=equationsToMatrix(k,q)

rank(S)

rank(C) %if rank(C)=rank([C,S]) the system can be solved

rank([C,S])

ji=C\S %ji is pi^(0) which i need to find and must fulfill all value less than 1 and non negative.

But, after i run, rank(C)/=rank([C,S]) so the system can't be solved

Is there any code that is wrong? or there is code that more simple than this code?

Answer 1

[Torsten]

Why don't you use

pi0 = null((A0+R*B).');
if ~isempty(pi0)
  pi0 = pi0.'/(pi0.'*inv(I-R)*e)
end

[DoinK]

@Torsten It shows

pi0 =

3×0 empty double matrix

It means no result right, Sir?

[Torsten]

Yes, it means (A0+R*B).' has full rank and thus null((A0+R*B).') is empty.

Thus obviously the same result as with your code.

[DoinK]

Oh i see, your code is the same as my code but shortened.

Thankyou Sir.

So there is one question left, how to find pi^(0).

I'm sure that matrix R is right, so it must be the code that is wrong (?)

[Torsten]

pi^0 does not exist for the matrices supplied.

I thought you had understood this because of the empty nullspace of (A0+R*B).' (or in your code: because rank(C)<rank([C,S])).

[DoinK]

Oh i get it, so the code already good.

The problems are the matrix (A_0, B, A, C, or R).

Thankyou so much Sir.

[DoinK]

I change my matrix, but still cannot find the pi^(0).

It is continuous time markov chain (sum of every row equal to zero),

Q=[A_0 C 0 0 0 ...; B A C 0 0 ...; 0 B A C 0 ...; ...] and the matrix only on tridiagonal of matrix Q.

(the sum of upper diagonal) / (the sum of lower diagonal) <1.

example the second row 3/12 <1

example the third row 3/(7+5)<1

all have the same value that is 1/4 < 1

Here are my matrix and matrix R by successive substitution.

A_0=[-3 3 0;12 -15 3;7 5 -15]

B=[1 6 5;0 4 4;0 0 7]

A=[-15 3 0;4 -15 3;3 2 -15]

C=[0 0 0;0 0 0;3 0 0]

mm=3;

R=zeros(mm);

X=0;

Z=0;

Err=0;

for i=1:500;

fprintf('i=%g',i);

R=((-C)-((R)^2)*B)*A^-1;

num2str(R,'%7.3f');

Err=R-X;

if isequal(R,X)

break

end

era=abs(Err);

erb=abs(X);

mera=max(era);

merax=max(max(era));

merb=max(erb);

merbx=max(max(erb));

ert=merax/merbx

if ert<=10^(-5)

break %iterasi i= 58

end

X=R;

num2str(X,'%7.3f');

end

From the linear equation pi^(0) * (A_0+R*B)=0,

because pi^(0) is not zero, so (A_0+R*B) must be zero.

[Torsten]

For that pi^0*M = 0 for a matrix M and a nonzreo vector pi^0, M does not need to be zero. pi^0 has to lie in the nullspace of M.'.

[DoinK]

Thankyou for your explanation Sir, it really help me.

One more question, from my matrix A_0, B, A, C, and R

If I compute manual

[a b c] * (A_0+RB)=[0 0 0] , I get 3 linear equation name 1, 2, and 3.

[a b c]*(inv(I-R))*[1 1 1].'=1 , I get 1 linear equation name 4.

From my calculation, If I use eq. 1 2 3 and 4, no solution exist.

If I use eq. 1 2 and 3, a=0, b=0, c=0.

If I use eq. 4 and 2 more equation between (1,2,or 3) , there are solutions for a, b, c but the value is not fix. Have a litte different if I use another equation to compare.

Example if I use eq. 4 1 2 compare to eq. 4 2 3, there is a difference in the 6th number after comma.

[Torsten]

I don't understand why you experiment above with solutions to 3 equations out of the 4.

For that your problem has a solution, there must exist a vector different from the null vector that satisfies pi^0*(A_0+R*B) = 0 or (A_0+R*B).'*pi^0.' = 0. If such vector(s) exist, you'll get a basis of the vector space spanned by these vectors by the command null((A_0+R*B).'). If the result of this command is empty, your problem has no solution.

Equivalently you can check whether rank((A_0+R*B).') < 3. If this is the case, you'll also be able to solve your problem.