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
______________________________________________ [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
