LWR model (deterministic)

lwrmodel.R

#1D
simLWR <- function(r, m, u0,
                   sys_noise=0, 
                   obs_noise=0,
                   explicit=FALSE){
  n <- length(u0)
  if(explicit){
    A <- diag(1-r,n)
    A[cbind(2:n,1:(n-1))] <- r    
  }else{
    A <- diag(1+r,n)
    A[cbind(2:n,1:(n-1))] <- -r    
    A <- solve(A)
  }
  b <- matrix(NA,m,n)
  b0 <- u0
  b[1,] <- u0
  for(i in 2:m){
    b1 <- A %*% b0
    b1 <- b1 + rnorm(n,0,sys_noise)
    b0 <- b1
    b[i,] <- b0
  }
  colnames(b) <- seq(0,by=1,length.out=n)
  rownames(b) <- seq(0,by=1,length.out=m)
  b <- b + rnorm(n*m,0,obs_noise)
  return(list(U=b, A=A))
}
#2D
simLWR_2d <- function(r, u0, timestep){
  m <- nrow(u0)
  n <- ncol(u0)
  A <- diag(-1,m)
  A[cbind(2:m,1:(m-1))] <- 1
  B <- diag(-1,n)
  B[cbind(1:(n-1),2:n)] <- 1
  U <- array(0, dim = c(m,n,timestep+1))
  U[,,1] <- u0
  for(i in 1:timestep){
    u1 <- u0 + r * (A %*% u0 + u0%*% B)
    u0 <- u1
    U[,,i+1] <- u1
  }
  return(list(U=U, A=A, B=B))
}

u0 <- matrix(0,5,5)
u0[1] <- 1

r <- 0.1
out <- simLWR_2d(r,u0,24)
dfU <- reshape2::melt(out$U)
colnames(dfU) <- c("x","y","time","density")
dfU$time <- dfU$time - 1L
####
#plot
library(ggplot2)
library(gganimate)
ggplot(dfU,aes(x=x,y=y,fill=density))+
  geom_tile()+
  facet_wrap(~time, labeller = label_both)+
  scale_fill_gradient(low="white", high="royalblue")+
  theme_classic()

ggsave("stil.png")
ggplot(dfU,aes(x=x, y=y, fill=density))+
  geom_tile()+
  scale_fill_gradient(low="white", high="royalblue")+
  theme_classic()+
  transition_time(time)+
  labs(title = 'time: {frame_time}')

anim_save("lwr.gif")