Re: MatrixPower[]

*To*: mathgroup at smc.vnet.net*Subject*: [mg28932] Re: [mg28919] MatrixPower[]*From*: andrzej <andrzej at bekkoame.ne.jp>*Date*: Sat, 19 May 2001 22:27:50 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

Undoubtedly yours it the best universal method and I hope something like it will be implemented in the next version. But I would also like to point out that in one can usually solve this problem in an "ad hoc" way as follows. First, perturb your matrix so it that it no longer has a zero eigenvalue, e.g., in your case you can take: M1 = {{1 + t, 1, 1}, {0, 2, 0}, {1, 1, 1}}; Now In[2]:= Simplify[Limit[MatrixPower[M1, n] + O[t], t -> 0], n > 0] Out[2]= {{2^(-1 + n), 2^(-1 + n)*n, 2^(-1 + n)}, {0, 2^n, 0}, {2^(-1 + n), 2^(-1 + n)*n, 2^(-1 + n)}} Which is the same answer as given by your method. Andrzej Kozlowski Toyama International University JAPAN http://platon.c.u-tokyo.ac.jp/andrzej/ http://sigma.tuins.ac.jp/~andrzej/ On Friday, May 18, 2001, at 02:13 PM, Konstantin L Kouptsov wrote: > > 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. > >