Home
PDF version of this note
Mathematica notebook

Symbolic solution to expansion of natural frequencies method for determining (e^At)

Written by Nasser M. Abbasi, based on lecture given by Professor Barmish, ECE 717, University of Wisconsin, Madison
11/1/2014

The method of finding (e^At) using expansion of natural frequencies is first applied to a 2×2 matrix to illustrate the symbolic method. Then a general purpose function is written at the end using the symbolic method and applied to much larger problems for illustrations. This is a standalone function which can be called to determine (e^At) for matrices A which can have repeated eigenvalues.  The result is verified by comparing the output to Mathematica buildin function MatrixExp[]

Example 1 from lecture done using symbolics

Define A

Clear["Global`*"]
A0={{2,0},{-3,-3}};
MatrixForm[A0]

Find eigenvalues of A

λ={2,-3};

Set the eigenvalue multiplers

multipliers={1,1}; (*multiplers of eigenvalues*)

Apply (e^At) algorithm using expansion of natural frequencies

eAt=Total@Flatten@Last@Reap@Do[Do[  Sow[Y[k,i] t^k Exp[t λ[[i]]]],{k,0,multipliers[[i]]-1}],{i,1,Length[λ]}];
Row[{"(e^At)=",eAt}]

Automatic generation of  lists of unknowns Y(k,i) and list of symbols to b[i] to use for solving Ax=b

varsb=Array[b,Total[multipliers],{0,Total[multipliers]-1}]
varsX=Flatten@Last@Reap@Do[Do[Sow@Y[k,i],{k,0,multipliers[[i]]-1}],{i,1,Length[multipliers]}]

Take repeated derivatives of (e^At) to generate the set of equations to solve for Y(k,i)

eqs=Table[varsb[[k+1]]==(D[eAt,{t,k}])/.t->0,{k,0,Total[multipliers]-1}]

The problem can be viewed as Ax=b where x=Y(k,i) is the unknowns and b are the symbolic names for the matrix powers.

{rightSide,A1}=Normal@CoefficientArrays[eqs,varsX];
Row[{MatrixForm[A1],MatrixForm[Transpose@{varsX}],"=",MatrixForm[rightSide]}]

Now solve Ax=b

sol=Solve[eqs,varsX]

Replace symbols b(i) on the right of the above, by the actual numerical values of the matrix powers

bN=(b[#]->MatrixPower[A0,#])&/@Range[0,Length[varsb]-1];
Row[{bN[[#,1]],"=",MatrixForm[bN[[#,-1]]]}]&/@Range[Length[bN]]

Now apply these numerical values to the solution found earlier in order to obtain numerical values for Y(k,i)

sol2=First[sol/.bN];
Row[{sol2[[#,1]],"=",MatrixForm@sol2[[#,-1]]}]&/@Range[Length[sol2]]

Now since Y(k,i) are now found,  evaluate (e^At) found above using the above solution

Row[{"(e^At)=",MatrixForm[(eAt/.sol2)]}]

Compare to Mathematica buildin function

MatrixExp[A0 t]//MatrixForm

Function to automatically calculate (e^At) using symbolic natural frequencies method

matrixExp[A0_,debug_,Y_,t_]:=Module[{eigenvalues,multipliers,tally,eAt,k,i,varsb,varsX,rightSide,A1,eqs,sol,sol2,bN},
eigenvalues=Eigenvalues[A0];
tally=Tally[eigenvalues];
eigenvalues=tally[[All,1]];
multipliers=tally[[All,2]];
If[debug,Print[Grid[Join[{{"λ","multiplier"}},tally],Frame->All]]];

eAt=Total@Flatten@Last@Reap@Do[Do[  Sow[Y[ k,eigenvalues[[i]] ] t^k Exp[t eigenvalues[[i]] ]],{k,0,multipliers[[i]] -1}],{i,1,Length[eigenvalues]}];

If[debug,Print@Row[{"(e^At)=",eAt}]];

varsb=Array[b,Total[multipliers],{0,Total[multipliers]-1}];
varsX=Flatten@Last@Reap@Do[Do[Sow@Y[k,eigenvalues[[i]]],{k,0,multipliers[[i]]-1}],{i,1,Length[multipliers]}];
eqs=Table[varsb[[k+1]]==(D[eAt,{t,k}])/.t->0,{k,0,Length[varsb]-1}];
If[debug,Print["eqs=\n",TableForm@eqs]];

{rightSide,A1}=Normal@CoefficientArrays[eqs,varsX];

If[debug,Print["A x = b  form"]];
If[debug,Print@Row[{MatrixForm[-A1],MatrixForm[Transpose@{varsX}],"=",MatrixForm[rightSide]}]];

sol=First@Solve[eqs,varsX];
bN=(varsb[[#+1]]->MatrixPower[A0,#])&/@Range[0,Length[varsb]-1];
sol2=sol/.bN;
If[debug,Print[Row[{sol2[[#,1]],"=",MatrixForm@sol2[[#,-1]]}]&/@Range[Length[sol2]]]];
eAt=eAt/.sol2;
If[debug,Print@Row[{"(e^At)=",MatrixForm[eAt,Dividers->All]}]];
eAt
]

Example 2 from lecture solved using the above function

Clear[Y,t];
A0={{1,0,0,0,0,0},{0,1,0,0,0,0},{-1,2,2,0,0,0},{1,0,-1,2,0,0},{0,0,1,0,2,0},{0,-1,1,0,0,3}};
(z1=matrixExp[A0,True,Y,t])//MatrixForm

Compare to buildin function to verify

(z2=MatrixExp[A0 t])//MatrixForm

MatrixForm[Simplify[z1-z2]]

Verification Method, suggested by Professor Barmish. We Differentiate Matrix (e^At) found above and see it is the same as (A^k) for t=0

tbl=Table[{k,MatrixForm[D[z1,{t,k}]/.t->0],MatrixForm[MatrixPower[A0,k]]},{k,0,5}];
h={{"k","((d^k)/dt)(e^At) at t=0","(A^k)"}};
Grid[Join[h,tbl],Frame->All]