The following is a quick example of a 'classical approach' To raising a
matrix to a power (integer or otherwise):
/* [wxMaxima: title start ]
Example of raising matrix to arbitrary power
[wxMaxima: title end ] */
/* [wxMaxima: comment start ]
First, set up the matrix...assume (2x2) Markov transition matrix
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
a : matrix([0.8,0.5],[0.2,0.5]);
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
Then, extract eigenvalues and eigenvectors....
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
[vals,vecs] : float(eigenvectors(a));
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
Transform eigenvectors from Maxima LISP-list to matrix form
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
e : transpose(apply('matrix,map('first,vecs)));
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
diagonalize matrix...
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
d : invert(e).a.e;
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
now, raise the matrix to the desire power - example is simple square
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
a2 : float(e.(d^2).invert(e));
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
now we raise it to a non-scalar power
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
a15 : float(e.(d^(1.5)).invert(e));
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
now we try square-root of matrix
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
roota : float(e.(d^(0.5)).invert(e));
/* [wxMaxima: input end ] */
/* [wxMaxima: comment start ]
The preceding is absolutely correct. The only major issue is whether or
not there is >1 square-root for a given matrix.
This is a non-trivial problem. I believe for non-negative
square-matrices, with all elements
0<=a<=1, I belive there is only a single square-root. Need to confirm
formally...
[wxMaxima: comment end ] */
/* [wxMaxima: comment start ]
Now we test, by taking the result and squaring it - should yield
original matrix
[wxMaxima: comment end ] */
/* [wxMaxima: input start ] */
float(roota.roota);
/* [wxMaxima: input end ] */