Re: e^A (*A is a matrix)
- To: mathgroup at yoda.physics.unc.edu
- Subject: Re: e^A (*A is a matrix)
- From: "Roger B. Kirchner" <kirchner at cs.umn.edu>
- Date: Tue, 17 Mar 92 07:05:46 -0600
MMa 2.0 has MatrixPower, but this can fail when a matrix doesn't have
a full set of eigenvectors. So too will any method that depends on
this.
The method in MatrixReplace.m just depends upon being able to find
the eigenvalues.
A couple of examples.
Mathematica 2.0 for NeXT
Copyright 1988-91 Wolfram Research, Inc.
-- NeXT graphics initialized --
In[1]:= << matrixreplace.m
In[2]:= A = {{2, 0}, {0, 3}}
Out[2]= {{2, 0}, {0, 3}}
In[3]:= MatrixReplace[E^(x t), x -> A]
2 t 3 t
Out[3]= {{E , 0}, {0, E }}
In[4]:= MatrixReplace[x^n, x -> A]
n n
Out[4]= {{2 , 0}, {0, 3 }}
The file:
(* MatrixReplace.m *)
(* Roger Kirchner, Carleton College, December 1991
rkirchne at mathcs.carleton.edu *)
(* MatrixReplace evaluates an expression at a matrix using a
generalization of a formula in
R. B. Kirchner, "An Explicit Formula for e^(A t)", American
Mathematical Monthly, Vol 74, No 10, December, 1967.
E.g., Let A be a square matrix.
MatrixReplace[f[x], x -> A] returns f[A], where multiplication
is matrix multiplication.
MatrixReplace[E^(x t), x -> A] returns the matrix exponential of A t.
MatrixReplace[x^n, x -> A] returns the nth matrix power of A.
Mathematica must be able to compute the eigenvalues of A.
*)
BeginPackage["MatrixReplace`"]
MatrixReplace::usage =
"MatrixReplace[expr, x -> mat] replaces x in expr by mat. Uses
eigenvalues of mat. Multiplication is matrix multiplication."
Begin["`Private`"]
taycoefs[f_, x_, x0_, n_] :=
Block[ {tp, u},
tp = Normal[Series[f, {x, x0, n}]];
CoefficientList[ tp /. x -> x0 + u, u]]
mathorner[coefs_, m_] :=
Block[{size = Length[m]},
If[coefs == {},
Table[0, {size}, {size}],
First[coefs] IdentityMatrix[size] +
m . mathorner[Rest[coefs], m]]]
zeros[poly_, x_] :=
x /. Solve[poly == 0, x]
duplist[lst_] :=
Block[{duplist1},
duplist1[ls_, x_, k_] :=
If[ls == {},
{{x, k}},
If[First[ls] === x,
duplist1[Rest[ls], x, k+1],
Prepend[duplist1[Rest[ls], First[ls], 1],
{x, k}]]];
duplist1[lst, First[lst], 0]]
zerolist[poly_, x_] :=
duplist[zeros[poly, x]]
van1[x_, s_, n_] :=
Block[{i, j},
Table[If[i < j,
Binomial[j-1, i-1] x^(j - i),
If[i == j,
1,
0]],
{i, 1, s}, {j, 1, n}]]
van2[pairs_, n_] :=
If[pairs == {},
{},
Block[{x, s},
x = pairs[[1, 1]];
s = pairs[[1, 2]];
Prepend[van2[Rest[pairs], n],
van1[x, s, n]]]]
vanmat[pairs_, n_] :=
Transpose[Flatten[van2[pairs, n], 1]]
rownums[cnts_] :=
Block[{sums},
sums[v0_, lst_] :=
If[lst == {},
{},
Prepend[sums[v0 + First[lst], Rest[lst]], v0]];
sums[1, cnts]]
qcoefs[pairs_, n_] :=
Block[{m, minv, cnts, rowstokeep},
m = vanmat[pairs, n];
minv = Inverse[m];
cnts = Map[#[[2]]&, pairs];
rowstokeep = rownums[cnts];
minv[[rowstokeep, Range[1, n]]]]
MatrixReplace[f_, x_ -> m_] :=
Block[{sumprods, eigenlist, eigens, qs, tays},
sumprods[qpart_, taypart_, eigens_] :=
Block[{size, id},
size = Length[m];
id = IdentityMatrix[size];
If[qpart == {},
Table[0, {size}, {size}],
mathorner[qpart[[1]], m] .
mathorner[taypart[[1]], m - eigens[[1]] id] +
sumprods[Rest[qpart], Rest[taypart],
Rest[eigens]]]];
eigenlist = duplist[Eigenvalues[m]];
eigens = Map[#[[1]]&, eigenlist];
qs = qcoefs[eigenlist, Length[m]];
tays =
Map[taycoefs[f, x, #[[1]], #[[2]] - 1]&,
eigenlist];
Simplify[sumprods[qs, tays, eigens]]]
End[]
Protect[ MatrixReplace ]
EndPackage[]
Roger Kirchner