Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am doing a work about Markov Chains and I have as task to determine, either by induction or by spectral methods, the form of the higher order transition matrix of the model.

My transition matrix (P) is:

transitionmatrix = {{17/72, 3/72, 5/72, 0, 0}, {5/72, 5/72, 7/72, 1/72, 0}, {4/72,9/72,
4/72, 0, 4/72}, {0, 0, 3/72, 1/72, 0}, {3/72, 1/72, 0, 0, 0}}      

By spectral method P=MDM^-1 :

In[6]:= eigenvalues = Eigenvalues[transitionmatrix] // N
Out[6]= {0.297046, 0.113193, -0.0215001 + 0.0136874 I, -0.0215001 - 0.0136874 I,    
0.00776109}

In[7]:= eigenvectors = Eigenvectors[transitionmatrix] // N
Out[7]= {{6.03283, 3.28881, 3.32027, 0.48858, 1.}, {4.22247, -4.51749, -4.76335, 
-1.99863, 1.}, {-0.171634 + 0.124637 I, -1.0331 + 0.61158 I, 1.23199 - 0.86313 I,     
-1.60369  + 0.395985 I, 1.}, {-0.171634 - 0.124637 I, -1.0331 - 0.61158 I, 1.23199 + 
0.86313 I, -1.60369 - 0.395985 I, 1.}, {0.280052, -0.281356, -0.752063, 5.11374, 1.}}

In[18]:= M = Transpose[eigenvectors]
Out[18]= {{6.03283, 
4.22247, -0.171634 + 0.124637 I, -0.171634 - 0.124637 I, 
0.280052}, {3.28881, -4.51749, -1.0331 + 0.61158 I, -1.0331 - 
0.61158 I, -0.281356}, {3.32027, -4.76335, 1.23199 - 0.86313 I, 
1.23199 + 0.86313 I, -0.752063}, {0.48858, -1.99863, -1.60369 + 
0.395985 I, -1.60369 - 0.395985 I, 5.11374}, {1., 1., 1., 1., 1.}}

In[16]:= Diag = DiagonalMatrix[eigenvalues]
Out[16]= {{0.297046 + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, 
0. + 0. I}, {0. + 0. I, 0.113193 + 0. I, 0. + 0. I, 0. + 0. I, 
0. + 0. I}, {0. + 0. I, 0. + 0. I, -0.0215001 + 0.0136874 I, 
0. + 0. I, 0. + 0. I}, {0. + 0. I, 0. + 0. I, 
0. + 0. I, -0.0215001 - 0.0136874 I, 0. + 0. I}, {0. + 0. I, 
0. + 0. I, 0. + 0. I, 0. + 0. I, 0.00776109 + 0. I}}

In[17]:= Minver = Inverse[M]
Out[17]= {{0.109069 + 7.70366*10^-19 I, 0.0487956 + 6.50521*10^-19 I, 
0.0514224 - 9.72945*10^-19 I, 0.00239343 + 5.5971*10^-19 I, 
0.00961738 - 3.82994*10^-18 I}, {0.0828826 - 
4.20783*10^-18 I, -0.0865795 + 1.0842*10^-17 I, -0.0549339 - 
8.27329*10^-18 I, -0.0121092 - 4.07692*10^-18 I, -0.0269617 + 
1.91893*10^-17 I}, {-0.0897774 + 0.141661 I, 
0.0551552 - 0.443078 I, -0.00936235 + 0.351851 I, -0.0773342 + 
0.143982 I, 0.429087 - 0.636007 I}, {-0.0897774 - 0.141661 I, 
0.0551552 + 0.443078 I, -0.00936235 - 0.351851 I, -0.0773342 - 
0.143982 I, 
0.429087 + 0.636007 I}, {-0.0123972 - 
4.92909*10^-18 I, -0.0725264 + 1.38778*10^-17 I, 
0.0222361 - 1.0567*10^-17 I, 0.164384 - 5.62511*10^-18 I, 
0.159171 + 2.19287*10^-17 I}}

But when i do

M.Diag.Minver

i dont get my transition matrix i get:

In[14]:= M.Diag.Minver

Out[14]= {{0.236111 - 7.80878*10^-19 I, 0.0416667 + 6.58726*10^-18 I, 
0.0694444 - 5.94759*10^-18 I, -3.43324*10^-18 - 8.88005*10^-19 I, 
4.15626*10^-17 + 2.63492*10^-18 I}, {0.0694444 + 2.15757*10^-18 I, 
0.0694444 - 3.80299*10^-18 I, 0.0972222 + 2.07205*10^-18 I, 
0.0138889 + 3.02264*10^-18 I, 
1.88651*10^-17 - 1.20859*10^-17 I}, {0.0555556 + 4.03989*10^-18 I, 
0.125 - 6.75901*10^-18 I, 
0.0555556 + 5.16039*10^-18 I, -1.15614*10^-17 + 2.29187*10^-18 I, 
0.0555556 - 1.62182*10^-17 I}, {-2.1617*10^-15 + 
2.2484*10^-20 I, -9.09499*10^-16 - 5.39346*10^-19 I, 
0.0416667 - 6.27859*10^-20 I, 
0.0138889 + 1.20315*10^-18 I, -2.11264*10^-16 - 
2.33486*10^-18 I}, {0.0416667 + 9.43055*10^-20 I, 
0.0138889 + 9.58975*10^-19 I, -1.44261*10^-17 - 
6.90605*10^-19 I, -2.78314*10^-18 - 5.28583*10^-19 I, 
1.27434*10^-17 + 4.46324*10^-19 I}}

Please could u help me, i dont know what i am doing wrong.

And could u help me finding the higher transition matriz by induction? My teacher gave this code, but i dont understand it, so i wouldnt like to use.

POTtrMx[0] = IdentityMatrix[Length[StSp]]; 
POTtrMx[n_Integer] := POTtrMx[n] = Dot[POTtrMx[n - 1], trMx]
share|improve this question

closed as off-topic by Szabolcs, Sjoerd C. de Vries, Öskå, rasher, RunnyKine Jun 10 at 0:02

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Szabolcs, Sjoerd C. de Vries, Öskå, rasher, RunnyKine
If this question can be reworded to fit the rules in the help center, please edit the question.

    
Matrix multiplication is done with ., not with *. If using the correct operator, you'll get the expected result. –  Szabolcs Jun 9 at 16:48
    
my mistake, wanted to put a dot. In mathematica I put a dot and i dind get the transition matrix. –  Mariana da Costa Jun 9 at 16:51
    
Please post a complete and reproducible example then, and show how you compared the matrices. I do get the same matrix, but your code is not complete: I had to define eigenvalues separate as its missing from your post. Maybe we did something differently. –  Szabolcs Jun 9 at 16:54
    
it is complete now. –  Mariana da Costa Jun 9 at 16:58
1  
OK, the output you posted is in fact the same as transitionmatrix to high precision. The difference between them is less than 10^-7 and it is due to numerical errors. The result you get is correct (up to a reasonable precision). –  Szabolcs Jun 9 at 17:41

1 Answer 1

Just to reiterate what was covered in the comments, any matrix $\mathbf{P}$ can be diagonalized via $$\mathbf{P}=\mathbf{M}\mathbf{D}\mathbf{M}^{-1}$$ where $$\mathbf{M}=\text{Eigenvectors}[\mathbf{P}]^\mathsf{T} \\\mathbf{D}=\text{Eigenvalues}[\mathbf{P}].$$

As you showed, this can be done by

P = N[{{17/72, 3/72, 5/72, 0, 0}, {5/72, 5/72, 7/72, 1/72, 0}, 
         {4/72, 9/72, 4/72, 0, 4/72}, {0, 0, 3/72, 1/72, 0}, 
         {3/72, 1/72, 0, 0, 0}}]; 
{d, Mt} = Eigensystem[P]; 
Max[Abs[Chop[
   P - Transpose[Mt] . DiagonalMatrix[d] . Inverse[Transpose[Mt]]]]]

Out: 0

which shows that they're equal.

To compute, for example, $\mathbf{P}^{5}$, there are three ways:

Transpose[Mt].DiagonalMatrix[d^5].Inverse[Transpose[Mt]] // Chop
MatrixPower[P, 5]
P.P.P.P.P

The method your teacher uses just computes it recursively, but it's sort of stupid, since there is already a built-in function to do that (MatrixPower).

share|improve this answer
    
Thanks for help :D –  Mariana da Costa Jun 9 at 20:31

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