I'm looking at this Mathematica tutorial where the author defines matrix function by summing over (matrix-valued) residues of the resolvent.
ClearAll["Global`*"];
d = 3;
A = {{1, 4, 16}, {18, 20, 4}, {-12, -14, -7}};
B = {{15, 35, 55}, {35, 35, 77}, {55, 77, 99}};
resolvent[s_, A_] := Inverse[s IdentityMatrix[d] - A];
resolventResidue[eig_, A_] :=
Simplify[(s - eig) resolvent[s, A]] /. s -> eig;
f[x_] = Exp[x];
MatrixFunction[f, A] ==
Total[f[#] resolventResidue[#, A] & /@ Eigenvalues[A]] (*True*)
MatrixFunction[f, B] ==
Total[f[#] resolventResidue[#, B] & /@ Eigenvalues[B]] (* False *)
However, his implementation relies on Simplify
which is not robust. IE, for matrix $B$, the expression is too complicated for Simplify to work. What is a better way to implement this?
resolventResidue[eig_, A_] := Limit[(s - eig) resolvent[s, A], s -> eig];
seems to work, applyN
to the result forB
matrix case $\endgroup$