Does anybody have any fixes for residue? I ask because ...
I have a function that computes the matrix exponential. I do this by
finding
the residues of the resolvent matrix at each of its poles (Laplace
inversion).
For simple tests, it works okay, but other tests reveal bugs in residue
and friends.
Here is an example
(C1) m : matrix([2,1,0],[1,2,1],[0,1,10])$
(C2) p : charpoly(m,z)$
Let's find the residues of 1/p at one of its poles. Each pole of 1/p is
simple:
(C3) resultant(p,diff(p,z));
(D3) 1
(C4) sp : solve(p,z)$
(C5) residue(1/p,z,rhs(sp[1]));
(D5) 0
This is bogus; this pole is simple, we should be able to find its residue
using the limit
(C6) limit((z-rhs(sp[1]))/p,z,rhs(sp[1]));
((42*SQRT(3)*SQRT(22)*%I+476)^(2/3)*((17*SQRT(3)*SQRT(22)-99
*SQRT(3))*%I+51*SQRT(22)+99)
+(525*SQRT(22)-175*SQRT(3)*SQRT(22)*%I)*(42*SQRT(3)*SQRT(22)*%I+476)^(1/3))
/646800
I believe this is correct. It seems that residue quickly trurns its work
over to taylor, so mostly
likely this is a taylor bug.
Note: my matrix exponential code (should) work for arbitrary square
matrices. It doesn't require that
the matrix be selfadjoint or normal, or ...
Barton