[R] matrix of eigenvalues

[Christian Jost]

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.

-- 
***********************************************************
http://cognition.ups-tlse.fr/vas-y.php?id=chj  jost at cict.fr
Christian Jost                                   (PhD, MdC)
Centre de Recherches sur la Cognition Animale
Universite Paul Sabatier, Bat IV R3
118 route de Narbonne
31062 Toulouse cedex 4, France
Tel: +33 5 61 55 64 37   Fax: +33 5 61 55 61 54