[R] matrix of eigenvalues
Christian Jost
jost at cict.fr
Wed Oct 20 16:26:11 CEST 2004
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
More information about the R-help
mailing list