Re:
- To: mathgroup at smc.vnet.net
- Subject: [mg28960] Re:
- From: "Allan Hayes" <hay at haystack.demon.co.uk>
- Date: Wed, 23 May 2001 01:54:28 -0400 (EDT)
- References: <okLN6.9484$dl4.29669@ralph.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Konstantin,
Here is some code that I have wrote several years ago, also using the
Cayley-Hamilton Theorem. I haven't had time to compare the code with your's,
but they must have a lot in common.
Numerical analysts have warned me that this is very tricky calculation, so I
would be interested to hear of any problems.
PolyMatrixReplace[poly_, var_, matr_] :=
With[{Id = IdentityMatrix[Length[matr]]},
(Fold[#1 . matr + #2*Id & , 0*matr, #1] & )[
Reverse[CoefficientList[poly, var]]]
]
MatrixReplace[expr_, x_ -> A_] :=
Module[{evals, order, eset, df, data, interpol},
evals = Eigenvalues[A]; order = Count[evals, #1] & ;
eset = Union[evals];
df = Derivative[#1][Function[x, expr]] & ;
data = ({#1, Table[df[k - 1][#1],
{k, order[#1]}]} & ) /@ eset;
interpol = InterpolatingPolynomial[data, x];
PolyMatrixReplace[interpol, x, A]
]
APPLICATION TO YOUR EXAMPLES
M1={{1,1,1},{0,1,0},{1,1,1}};
MatrixReplace[x^n,x\[Rule]M1]//Simplify
{{(1/2)*(0^n + 2^n), -1 + 2^n, (1/2)*(-0^n + 2^n)},
{0, 1, 0}, {(1/2)*(-0^n + 2^n), -1 + 2^n,
(1/2)*(0^n + 2^n)}}
%/. 0^n ->0
{{2^(-1 + n), -1 + 2^n, 2^(-1 + n)}, {0, 1, 0},
{2^(-1 + n), -1 + 2^n, 2^(-1 + n)}}
M2={{1,1,1},{0,2,0},{1,1,1}};
MatrixReplace[x^n, x->M2]//Simplify
{{(1/2)*(0^n + 2^n), 2^(-1 + n)*n, (1/2)*(-0^n + 2^n)},
{0, 2^n, 0}, {(1/2)*(-0^n + 2^n), 2^(-1 + n)*n,
(1/2)*(0^n + 2^n)}}
%/. 0^n->0
{{2^(-1 + n), 2^(-1 + n)*n, 2^(-1 + n)}, {0, 2^n, 0},
{2^(-1 + n), 2^(-1 + n)*n, 2^(-1 + n)}}
Allan
---------------------
Allan Hayes
Mathematica Training and Consulting
Leicester UK
www.haystack.demon.co.uk
hay at haystack.demon.co.uk
Voice: +44 (0)116 271 4198
Fax: +44 (0)870 164 0565
> 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.
>
>