Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1992
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1992

[Date Index] [Thread Index] [Author Index]

Search the Archive

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





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