gganimate; discretize heat equation

difusion2d.R

difueq <- function(u_ini, timestep){
  Deltat <-1/400
  Deltax <-1/10
  Deltay <-1/10
  r <- Deltat/(Deltax^2+Deltay^2)
  M <- nrow(u_ini)
  N <- ncol(u_ini)
  u_old <- u_ini
  U <- array(0, dim=c(M,N,timestep+1))
  U[,,1] <- u_ini
  tmp <- matrix(0,M,N)
  for(t in 1:timestep){
    for (i in 2:(M-1)) {
    for(j in 2:(N-1)){
      tmp[i,j] <- u_old[i,j+1] -  2*u_old[i,j] + u_old[i,j-1] +
        u_old[i+1,j] -  2*u_old[i,j] + u_old[i-1,j] 
    }
    }
    u_new <- u_old + r * tmp
    u_new[1,] <- u_new[2,]
    u_new[M,] <- u_new[M-1,]
    u_new[,1] <- u_new[,2]
    u_new[,N] <- u_new[,N-1]
    u_old <- u_new
    U[,,t+1] <- u_new
  }
  return(U)  
}

u_ini <- matrix(0,20,20)
u_ini[11:15,11:15] <- 1
u_ini[13,15] <- 0 
u_ini[14:15,11] <- 0 
u_ini[11:12,11] <- 0 
u_ini[15,12:13] <- 0 
u_ini[11,12:13] <- 0 
#image(u_ini)
U <- difueq(u_ini,50)


outdf <-reshape2::melt(U)
colnames(outdf) <- c("x","y","time","value")
head(outdf)

library(gganimate)
ggplot(outdf)+
  geom_tile(aes(x=x,y=y,fill=value))+
  scale_fill_gradient(low="white", high="royalblue")+
  transition_time(time)+
  theme_bw()

anim_save("difu.gif")