[Date Index]
[Thread Index]
[Author Index]
MatrixPower[]
*To*: mathgroup at smc.vnet.net
*Subject*: [mg28919] MatrixPower[]
*From*: Konstantin L Kouptsov <klk206 at nyu.edu>
*Date*: Fri, 18 May 2001 01:13:26 -0400 (EDT)
*Organization*: New York University
*Sender*: owner-wri-mathgroup at wolfram.com
Hello, Mathgroup,
I needed to calculate the n-th power of a given matrix M.
Usually it can be done with the MatrixPower[M,n] function.
But for the symbolic power it becomes a bit tricky.
Take M={{1,1,1},{0,1,0},{1,1,1}},
then
MatrixPower[M,n] gives
{{2^(-1+n),-1+2^n,2^(-1+n)},
{0, 1,0},
{2^(-1+ n),-1+2^n,2^(-1+n)}}
Now, if you change the matrix to
M={{1,1,1},{0,2,0},{1,1,1}}, then the MatrixPower function
gives an error
MatrixPower::"zvec": "Cannot compute MatrixPower because the
argument has a zero eigenvector."
which seems irrelevant, since the power of a square matrix always
exist.
The reason for this error is in the way Mathematica calculates the
power. Here is the explicit form of the used algorithm
{eval, evec} = Eigensystem[M0];
s = Transpose[evec];
A = DiagonalMatrix[eval];
M[n_] := s.MatrixPower[A, n].Inverse[s]
which of course does not work if s has no inverse.
To circumvent the problem one can use the Hamilton-Cailey theorem
(which can be found in the Korn,Korn mathematical handbook).
The theorem gives a recipe to calculate any function of a matrix
using Vandermonde matrix built on its eigenvalues and a finite
number of powers of the matrix itself.
Below is the implementation of this theorem. There is one thing to watch
out: it the matrix has equal eigenvalues, the Vandermonde matrix is
zero, so you get indeterminate 0/0. The way the theorem can be applied
is that you can resolve indeterminancy
by cancelling out zeroes in numerator and denominator, as in
(2-2)(7-5) 1
__________ = _
(2-2)(9-3) 3
which must be done in mathematically correct way (remember calculus?)
I used the function ReplaceEqual[] to discern the similar eigevalues
in the following manner:
ReplaceEqual[{1,1,2,3,4,4,5}] gives
{{1,1+i1,2,3,4,4+i2,5},{i1,i2}}
i.e. shift similar eigenvalues to make them different, and at the end
calculate the result by setting i1,i2 to zero. This is done by the following
code
ReplaceEqual[x_List] := Module[{i = 1, a, k, ls1 = {}, ls2, s, x1},
x1 = Sort[x]; a = First[x1]; ls2 = List[a];
For[k = 2, k <= Length[x1], k++,
If[x1[[k]] == a,
s = ToExpression["i" <> ToString[i++]];
AppendTo[ls2, x1[[k]] + s]; AppendTo[ls1, s],
AppendTo[ls2, x1[[k]]]; a = x1[[k]]
]
];
{ls2, ls1}
]
Now the Vandermonde matrix determinant can be calculated in a rather
efficient way
VandermondeMatrix[eval_List] := Module[{x, y},
x = Table[1, {Length[eval]}];
y = (eval^# &) /@ Range[1, Length[eval] - 1];
PrependTo[y, x]
]
but we also need the determinant of the Vandermonde matrix with one of
the columns replaced by {F[eval[[1]]],F[eval[[2]]],...}
VandermondeMatrixReplaced[eval_List,f_, n_] := Module[{w, fs},
w = VandermondeMatrix[eval]; (* or pull this call out for efficiency *)
fs = Map[f, eval];
ReplacePart[w, fs, n]
]
Calculation of this deteminant was a problem which caused my recent messages to
MathGroup. Rasmus Debitsch have suggested using the row echelon form
of the matrix. For smaller matrices simple Det[] works.
Finally we arrive at the HamiltonCailey[] function which computes any
function of a square matrix:
HamiltonCailey[A_, F_] := Module[{eig, n, w, wr, x},
{eig, vars} = ReplaceEqual[Eigenvalues[A]];
n = Dim[A]; (* dimension of the n x n matrix A *)
w = VandermondeMatrix[eig];
wr = (1/Det[w]) Sum[
Det[VandermondeMatrixReplaced[eig, F, n - k + 1]]*
MatrixPower[A, n - k], {k, 1, n}]
]
Do not forget to apply the Limit[] to the result. There is still a lot
to be said about the above method, but you get the idea.
For the matrix in above example
M={{1,1,1},{0,2,0},{1,1,1}}
you write
Limit[HamiltonCailey[M,#^n&]//Simplify,i1->0]
to get the answer
{{2^(-1+n),n 2^(-1+n),2^(-1+n)},
{0,2^n,0},
{2^(-1+n),n 2^(-1+n),2^(-1+n)}}
Konstantin.
PS. I wish to thank Rasmus Debitsch and all who responded to my post.
Prev by Date:
**Re: Creating graph with only a few data points**
Next by Date:
**Request for contour plotting algorithm**
Previous by thread:
**Re: Creating graph with only a few data points**
Next by thread:
**Re: MatrixPower[]**
| |