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 1eigen(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 [EMAIL PROTECTED] 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
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
