Numerical solution of the diffusion equation (finite difference approximation )

diff2d.r

library(gganimate)
library(ggplot2)
#install.packages("gganimate")
diffeq <- function(u_ini, r, timestep){
  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[9:13,9:13] <- 1
u_ini[11,13] <- 0 
u_ini[12:13,9] <- 0 
u_ini[9:10,9] <- 0 
u_ini[13,10:11] <- 0 
u_ini[9,10:11] <- 0 
#image(u_ini)

U <- diffeq(u_ini,0.2,55)

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

ggplot(outdf)+
  geom_tile(aes(x=x,y=y,fill=value))+
  scale_fill_gradient(low="white", high="royalblue")+
  facet_wrap(~time)+
  theme_classic()
ggsave("still.png")

head(outdf)
ggplot(outdf,aes(x=x,y=y,fill=value))+
  geom_tile()+
  scale_fill_gradient(low="white", high="royalblue")+
  theme_classic()+
  transition_time(time)+
  labs(title = 'step: {frame_time}')
anim_save("difu.gif")