swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: a9b4d5f19cb88cd9d0d31fdde0d654737327d085 authored by HwB on 15 April 2012, 00:00:00 UTC
version 1.0.5
version 1.0.5
Tip revision: a9b4d5f
expm.R
##
## e x p m . R Matrix Exponential
##
expm <- function(A, np = 128) {
if (!is.numeric(A) || !is.matrix(A) || nrow(A) != ncol(A))
stop("Argument 'A' must be a square numeric matrix.")
if (!is.numeric(np) || length(np) != 1 ||
floor(np) != ceiling(np) || np < 2)
stop("Argument 'np' must be an integer greater or equal to 2.")
N <- nrow(A)
circle <- exp(2i*pi*(1:np)/np) # generate np unit roots
z0 <- ceiling(mean(range(real(eig(A)))) + 0.1)
radius <- ceiling(max(abs(eig(A) - z0)) + 0.1)
z <- z0 + radius*circle
I <- eye(N); B <- zeros(N)
for (i in 1:np) {
R <- inv(z[i]*I - A) # resolvent matrix at point z(i)
B <- B + R * (z[i]-z0) * exp(z[i]) # add up contributions to integral
}
B <- real(B)/np
return(zapsmall(B))
}
logm <- function(A) {
if (!is.numeric(A) || !is.matrix(A) || nrow(A) != ncol(A))
stop("Argument 'A' must be a square numeric matrix.")
E <- eigen(A)
e <- E$values
if (any(imag(e) == 0 && real(e) <= 0))
stop("A must not have any nonpositive real eigenvalues.")
D <- diag(log(E$values))
X <- E$vectors %*% D %*% solve(E$vectors)
return(real(X))
}