Displaying 1 result from an estimated 1 matches for "res_eigen".
2010 Jul 17
0
Computing the power of a matrix efficiently and accurately
...(500,20,2)) # set up tridiagonal matrix
for(i in 1:(m-1))
{
X[i,i+1]=rnorm(1,5,1)
X[i+1,i]=rnorm(1,5,1)
}
n=100 # the power of the matrix X i.e. X^n = X%*%X%*%X ... n times (n = 1000
in my application)
X_eigen=eigen(X)
P=X_eigen$vectors
P_inv=solve(P)
d=X_eigen$values
M=diag(m)
res_mult=c()
res_eigen=c()
for(i in 1:n)
{
# straightforward multiplication...
M = M%*%X
res_mult = rbind(res_mult, M[m,])
# eigen-decomposition...
res_eigen = rbind(res_eigen, (P%*%diag(d^i)%*%P_inv)[m,])
}
# Look how the diagonal element X[m,m] changes with successive powers:
> res_mult[1,m]
[1] 23.45616
> res_...