deSolve ode.2D in R - changing parameter in spatial dimensions

Question

I am trying to solve a model of ordinary differential equations in 2D using the function ode.2D in the package deSolve in R. My current code is:

library(deSolve)
ModelDif <- function (time, y, pars, N, Da, dx) {
NN <- N*N
Nsp <- matrix(nrow = N, ncol = N,y[1:NN])

with (as.list(pars), {
dNsp <- Rmax * Nsp * (1- Nsp/K)
         zero <- rep(0, N)
         FluxNsp <- -Da * rbind(zero,(Nsp[2:N,] - Nsp[1:(N-1),]), zero)/dx
         dNsp    <- dNsp - (FluxNsp[2:(N+1),] - FluxNsp[1:N,])/dx
         FluxNsp <- -Da * cbind(zero,(Nsp[,2:N] - Nsp[,1:(N-1)]), zero)/dx
         dNsp    <- dNsp - (FluxNsp[,2:(N+1)] - FluxNsp[,1:N])/dx
         return(list(as.vector(dNsp)))
      })
     }
        
pars    <- c(Rmax  = 1.0,K=5)  
     R  <- 20                  
     N  <- 50                  
     dx <- R/N                   
     Da <- 0.03           
     NN <- N*N                    
     yini    <- rep(0, N*N)
     cc      <- c((NN/2):(NN/2+1)+N/2, (NN/2):(NN/2+1)-N/2)
     yini[cc] <- 1
     
     tiempo   <- seq(0, 20, by = 1)
     Final <- ode.2D(y = yini, times = tiempo, func = ModelDif, parms = pars,dimens = c(N, N), names = c("Nsp"),N = N, dx = dx, Da = Da, method = rkMethod("rk45ck"))

##Plotting
     
P1 <- subset(Final, select = "Nsp", arr = TRUE)
for(i in 1:length(tiempo)){
TT<-tiempo[i]
image(as.matrix(P1[,,i]), xlab = "Lat", ylab = "Lon")
mtext(paste("tiempo = ", TT),side=3)
Sys.sleep(0.5)}

However, now I need to replace the fixed parameter K in the main equation by the value obtained from the respective cell in a matrix with the same size of the grid used for solving the model.

For example, replace the fixed value of K by K<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N) and at each evaluation of the model replace K in the equation by the value in the corresponded cell.

I have search for solution in the site, but none has worked. I'll appreciate any advice. Thanks

score 1 · Accepted Answer · answered Oct 02 '20 at 18:26

If you use K<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N), then you need to replace all zeros with a positive number as you are using it as a denominator.

In the model below, I am replacing 0 with 1. You can replace it with another number of your choice.

ModelDif <- function (time, y, pars, N, Da, dx) {
  NN <- N*N
  Nsp <- matrix(nrow = N, ncol = N,y[1:NN])
  KK<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N)
  for (i in 1:NN) {
    if (round(KK[i],1) == 0) {KK[i] <- 1}  ###  replace all zeros with 1
  }
  NspK <- Nsp/KK
  
  with (as.list(pars), {
    dNsp <- Rmax * Nsp * (1- NspK)
    zero <- rep(0, N)
    FluxNsp <- -Da * rbind(zero,(Nsp[2:N,] - Nsp[1:(N-1),]), zero)/dx
    dNsp    <- dNsp - (FluxNsp[2:(N+1),] - FluxNsp[1:N,])/dx
    FluxNsp <- -Da * cbind(zero,(Nsp[,2:N] - Nsp[,1:(N-1)]), zero)/dx
    dNsp    <- dNsp - (FluxNsp[,2:(N+1)] - FluxNsp[,1:N])/dx
    return(list(as.vector(dNsp)))
  })
}

Also, it takes a while to run and gives the following warning: Warning in rk(y, times, func, parms, method = method, ...) : Number of time steps 104651 exceeded maxsteps at t = 16.6229

You may need to further adjust some parameters.

deSolve ode.2D in R - changing parameter in spatial dimensions

1 Answers1