Re: [R] matrix of eigenvalues

[Spencer Graves]

     That sounds like what the default method of "mexp" does, though I 
haven't checked it to be sure.  That works fine if the eigenvalue part 
of the Jordan canonical form is diagonal.  That does not hold for your 
examples.  However, the theorem I mentioned from one of Bellman's books 
says that any matrix can be approximated arbitrarily closely with 
another matrix with unique eigenvalues, for which the default method of 
"mexp" seems to work.  For more, see the classic paper by Moler and van 
Loan that Doug mentioned on "Nineteen Dubious Ways to Calculate the 
Matrix Exponential". 

     hope this helps. 
     spencer graves

Christian Jost wrote:
At 12:12 -0500 21/10/04, Douglas Bates wrote:
Spencer Graves wrote:
     I think you need the Schur decomposition, which seems currently 
not to be available in R.  The documentation for the the Matrix 
package describes a "Schur" function, but it's not available in the 
current Matrix package, as mentioned in my post on this issue 
yesterday (subj:  Schur decomposition).
     Moreover, Lindsey's mexp in rmutils won't work either:
 > A = matrix(cbind(c(-1,1),c(-4,3)),nrow=2)
 > mexp(A)
Error in solve.default(z$vectors) : system is computationally 
singular: reciprocal condition number = 4.13756e-017
 > mexp(A, "series")
Error in t * x : non-numeric argument to binary operator
 >
     Have you considered adding a little noise: > mexp(A+1e-6*rnorm(4))
         [,1]       [,2]
[1,] -2.718284 -10.873139
[2,]  2.718284   8.154859
 > mexp(A+1e-6*rnorm(4))
                       [,1]                     [,2]
[1,] -2.718284-1.060041e-12i -10.873126-1.575184e-12i
[2,]  2.718284+7.492895e-13i   8.154846+1.578515e-12i
 > mexp(A+1e-6*rnorm(4))
                       [,1]                     [,2]
[1,] -2.718284+6.146195e-13i -10.873127+1.586731e-12i
[2,]  2.718284-4.004574e-13i   8.154847-8.515411e-13i
 > mexp(A+1e-6*rnorm(4))
                       [,1]                     [,2]
[1,] -2.718283+3.077538e-13i -10.873130+2.433609e-13i
[2,]  2.718283-3.197442e-14i   8.154847-2.782219e-13i
 >
     hope this helps.  spencer graves
(I believe this is described in one of Richard Bellman's matrix 
analysis book.)

I rather suspected that that original question was about the matrix 
exponential and linear systems of differential equations.  The 
solution to such a system can only be written using the matrix 
exponential for diagonalizable systems and this is the classic 
example of a nondiagonalizable system.

Calculation of the matrix exponential is an operation that seems 
straightforward in theory and can be very difficult in practice. 
Moler and van Loan have a classic paper on "Nineteen Dubious Ways to 
Calculate the Matrix Exponential" which I would recommend reading.

interesting, I would not have suspected such difficulties based on the 
technical definition of the matrix exponential I know
exp(A) = Identity + A + A^2/2! + A^3/3! ...
which does not at all require A to be diagonalizable.

Now, the original question was indeed how to solve numerically a 
linear ODE system with non-diagonalizable matrix. The example I gave 
in another reply suggests that this could be done by taking all 
linearly independent eigenvectors, complement them to a full base P, 
then compute D=invP %*% A %*% P (upper triangular matrix), decompose 
it into a diagonal matrix + an upper triangular matrix with 0 on the 
diagonal for which exp(D) would be rather straightforward to compute 
as long as the system is not too large.

What I do not know is whether the complementing to a full base is 
feasible. Any idea?

Thanks, Christian.
--
Spencer Graves, PhD, Senior Development Engineer
O:  (408)938-4420;  mobile:  (408)655-4567
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html