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
-- 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

```

• Prev by Date: Interesting notebooks?
• Next by Date: Re: A definite Integral
• Previous by thread: Re: e^A (*A is a matrix)