[R] matrix of eigenvalues

[Spencer Graves]

      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.) 

Christian Jost wrote:

> Thanks for all the replys. I understand now that eigen does not return 
> a matrix with linearly independent eigenvectors (which we would need 
> to get its inverse).
>
> The question whether there is a function in R that does so remains. As 
> an example, consider the problem
> dx1/dt = -x1 - 4 x2
> dx2/dt =  x1 + 3 x2
> corresponding to the matrix system
> dx/dt = A x, x(0) = x0 = c(1,2) (for example)
>
> Analysing this system in R gives
> A = matrix(cbind(c(-1,1),c(-4,3)),nrow=2)
> for which eigen() returns a singular eigenvector matrix
> det(eigen(A)$vectors)
>
> Now, if I tried to solve the system by hand, I would find one 
> eigenvalue of value 1 and an associated eigenvector v1=[2, -1]. I can 
> complement this eigenvector to get a base, e.g. v2=[0,1] (which is 
> linearly independent of v1). Defining now a passage matrix as
> P = matrix(cbind(c(2,-1),c(0,1)),nrow=2)
> I can compute the inverse
> invP = solve(P)
> and
> invP %*% A %*% P -> D
> will return an upper triangular matrix (with eigenvalues on the 
> diagonal) for which the exponential can easily be computed by hand
>
> expDt = diag(exp(t),2,2) %*% matrix(c(1,0,-2*t,1),nrow=2)
> and the solution of the original system becomes
>
> x(t) = P %*% expDt %*% invP %*% x0
> and I can simulate a trajectory
>
> times = seq(0,2,by=0.1);
> x0 = c(1,2);
> sol = matrix(0,nrow=2,ncol=length(times));
> for (i in 1:length(times)) {
>   t = times[i];
>   expDt = diag(exp(t),2,2) %*% matrix(c(1,0,-2*t,1),nrow=2);
>   sol[,i] = P %*% expDt %*% invP %*% x0
> }
> plot(times,sol[1,],type='l',ylim=c(min(sol),max(sol)))
> lines(times,sol[2,],lty='dashed')
>
> Thanks to all who followed down here. Now, is there a function in R 
> that will construct my P automatically? Or should I apply directly the 
> function mexp(A) (from package mrutil) and completely bypass the 
> passage matrix stuff (which my math teachers like so much ;-)
>
> Best, Christian.
>
>> Christian Jost wrote:
>>
>>> I thought that the function
>>> eigen(A)
>>> will return a matrix with eigenvectors that are independent of each 
>>> other (thus forming a base and the matrix being invertible). This 
>>> seems not to be the case in the following example
>>> A=matrix(c(1,2,0,1),nrow=2,byrow=T)
>>> eigen(A) ->ev
>>> solve(ev$vectors)
>>>
>>> note that I try to get the upper triangular form with eigenvalues on 
>>> the diagonal and (possibly) 1 just atop the eigenvalues to be used 
>>> to solve a linear differential equation
>>> x' = Ax, x(0)=x0
>>> x(t) = P exp(D t) P^-1 x0
>>> where D is this upper triangular form and P is the "passage matrix" 
>>> (not sure about the correct english name) given by a base of 
>>> eigenvectors. So the test would be
>>> solve(ev$vectors) %*% A %*% ev$vectors - D
>>> should be 0
>>>
>>> Thanks for any help, Christian.
>>>
>>> ps: please copy reply also to my address, my subscription to the 
>>> R-help list seems to have delays
>>
>>
>> That particular matrix has repeated eigenvalues and a degenerate 
>> eigenspace.
>>
>>>  A <- matrix(c(1,0,2,1),nc=2)
>>>  A
>>
>>      [,1] [,2]
>> [1,]    1    2
>> [2,]    0    1
>>
>>>  eigen(A)
>>
>> $values
>> [1] 1 1
>>
>> $vectors
>>      [,1]          [,2]
>> [1,]    1 -1.000000e+00
>> [2,]    0  1.110223e-16
>
>
>

-- 
Spencer Graves, PhD, Senior Development Engineer
O:  (408)938-4420;  mobile:  (408)655-4567