[R] matrix of eigenvalues

[Christian Jost]

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


-- 
***********************************************************
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