[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