So I'm in a stochastic simulations class and I having issues with the amount
of time it takes to run the Ising model.
I usually don't like to attach the code I'm running, since it will probably
make me look like a fool, but I figure its the best way I can find any bits
I can speed up run time.
As for the goals of the exercise:
I need the state of the system at time=1, 10k, 100k, 1mill, and 10mill
and the percentage of vertices with positive spin at all t
Just to be clear, i'm not expecting anyone to tell me how to program this
model, cause I know what I have works for this exercise, but it takes far
too long to run and I'd like to speed it up by replacing slow operations
wherever possible.
Thanks in advance for any help.
Code Follows:
N <- 200 # Size of model
R <- 10^7 # Runs
T <- 1 # Temperature
S <- matrix( sample(c(-1,1),N^2,replace=TRUE) ,nrow=N,ncol=N)
M <- cbind(rep(0,N+2),rbind( rep(0,N), S, rep(0,N)), rep(0,N+2))
q <- 1
Out <- c()
Out[[q]] <- M
pos <- rep(0,R)
pos[1] <- sum(M==1)
for(k in 1:R){
## Pick random vertex
U <- sample(0:(N^2-1),1)
i <- floor(U/N) + 2
j <- U%%(N) + 2
## Calculate Energy
Nei <- c(M[i-1,j],M[i+1,j],M[i,j-1],M[i,j+1])
Ei <- 2 * sum( M[i,j] != Nei)
Ej <- 2 * sum( sign(M[i,j])*(-1) != Nei)
## Accept Criteria Check
if( Ej < Ei ){
M[i,j] <- sign(M[i,j])*(-1)
pos[k+1] <- M[i,j]
} else {
if( runif(1) < exp(-1/T*(Ej-Ei)) ){
M[i,j] <- sign(M[i,j])*(-1)
pos[k+1] <- M[i,j]
} }
## Store state at time 10^4, 10^5, 10^6, 10^7
if( k %in% c(10^4,10^5,10^6,10^7) ){
q <- q+1
Out[[q]] <- M
}
}
## Output
image(Out[[1]])
image(Out[[2]])
image(Out[[3]])
image(Out[[4]])
image(Out[[5]])
plot( cumsum(pos)/N^2, pch='.')
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
