H i, all!
First of all, I'd like to apologize for my poor English. It's for years
I don't use it.
This is a R-version of a function I wrote a long ago for my HP48 calculator.
It works with the binary expression of the power and just need to
duplicate the mem used by X.
Hope this helps.
mtx.exp<-function(X,n)
#Function to calculate the n-th power of a matrix X;
{
phi <- diag(rep(1,length(X[1,])))
pot <- X #This is the first power of the matrix.
while (n > 0)
{
if (n%%2)
{
phi <- phi%*%pot;
}
n <- n%/%2;
pot <- pot %*% pot;
}
return(phi);
}
#Here is some output:
> xx <- matrix(c(1,0,1,1),2,2)
> xx
[,1] [,2]
[1,] 1 1
[2,] 0 1
> mtx.exp(xx,3)
[,1] [,2]
[1,] 1 3
[2,] 0 1
> mtx.exp(xx,4)
[,1] [,2]
[1,] 1 4
[2,] 0 1
> mtx.exp(xx,6)
[,1] [,2]
[1,] 1 6
[2,] 0 1
> mtx.exp(xx,10)
[,1] [,2]
[1,] 1 10
[2,] 0 1
> mtx.exp(xx,1000)
[,1] [,2]
[1,] 1 1000
[2,] 0 1
Vicente D. Canto Casasola
>>>>> "Vicente" == Vicente Canto Casasola <vicented.canto.ext at juntadeandalucia.es> >>>>> on Thu, 22 Jan 2004 14:34:01 +0100 writes:Vicente> H i, all! First of all, I'd like to apologize for Vicente> my poor English. It's for years I don't use it. no problem at all. Vicente> This is a R-version of a function I wrote a long Vicente> ago for my HP48 calculator. It works with the Vicente> binary expression of the power and just need to Vicente> duplicate the mem used by X. excellent. This is really the way to solve the problem I think. As I've mentioned earlier in this thread, computing a matrix "power" is really much easier than the matrix exponential. Hence I wouldn't use exponential in the function name. Also note that trailing ";" are considered as `dirty' (they are completely superfluous). These slight modifications (+ initial "test") give matPower <- function(X,n) ## Function to calculate the n-th power of a matrix X { if(n != round(n)) { n <- round(n) warning("rounding exponent `n' to", n) } phi <- diag(nrow = nrow(X)) pot <- X # the first power of the matrix. while (n > 0) { if (n %% 2) phi <- phi %*% pot n <- n %/% 2 pot <- pot %*% pot } return(phi) } Regards, Martin Maechler <maechler at stat.math.ethz.ch> http://stat.ethz.ch/~maechler/ Seminar fuer Statistik, ETH-Zentrum LEO C16 Leonhardstr. 27 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-1-632-3408 fax: ...-1228 <><