[R-sig-dyn-mod] Adding a delay to a ode.2d model

[Raffaello Vardavas]

Hello. I wonder if anyone can help me with the following.  I've modified the example code in the documentation for the function tran.2D so that I have varying drift and diffusion terms. This works and my code below runs without problems. Its just an example model - I'm in the exploration phase of the possibilities.


Now I'm attempting to put a lag in the my PDE. I understand that lags can be implemented as long as the solver is specified to be one of those comparable with time lags as given in the table found in the slides http://desolve.r-forge.r-project.org/slides/tutorial.pdf. I'm specifying the solver to be lsodes - which is one of the compatible solvers. Nevertheless, I'm getting problems when I begin to explore implementing lag terms in the Diff2Dc function shown in the code below. In particular, when I un-comment the part of the code that reads


  # if(t-delay>0){
  #   #tmp<-lagvalue(t-delay)
  #   browser()
  # }


I get an error message. Basically the function lagvalue does seem to work in this implementation - despite the fact I'm actually not using it for the time being to affect the PDE. Other examples of using the lagvalue when integrating a 1D model seems to work.


Please let me know what is wrong and how does one go about introducing a lag in a 2D PDE using the method of lines.


Thank you.

Raff



Code:


library(ReacTran)

#TODO: check "Fitting ODE in R, with the use of the FME package"

N     <- 51          # number of grid cells
ini   <- 1           # initial value at x=0
N2    <- ceiling(N/2)

dx <- 10/N
dy <- 10/N


### Diffusion
gloD.grid <- list()
# Diffusion on x-interfaces
gloD.grid$x.int <- matrix(nrow = N+1, ncol = N,
                          data = 2*runif(N*(N+1)))
# Diffusion on y-interfaces
gloD.grid$y.int <- matrix(nrow = N, ncol = N+1,
                          data = 1.5*runif(N*(N+1)))


### Drift
gloDrift.grid <- list()
# Diffusion on x-interfaces
gloDrift.grid$x.int <- matrix(nrow = N+1, ncol = N,
                          data = -0.3)
# Diffusion on y-interfaces
gloDrift.grid$y.int <- matrix(nrow = N, ncol = N+1,
                          data = -0.1)

# birth rate
inflow.grid <- matrix(nrow = N, ncol = N, data = 0)
inflow.grid[40, 40] <- 0.01

# outflow rate
outflow.rate.grid <- matrix(nrow = N, ncol = N, data = 0.1/(N*N))


### Parms
parms<-list(inflow.grid =inflow.grid,
            outflow.rate.grid=outflow.rate.grid,
            v.grid=gloDrift.grid,
            D.grid=gloD.grid)


# The model equations
Diff2Dc <- function (t, y, parms,CONCini,delay)  {
  with(as.list(parms),{

  CONC  <- matrix(nrow = N, ncol = N, data = y)

  dCONC <- tran.2D(CONC, dx = dx, dy = dy,
                   v.grid=v.grid,
                   D.grid = D.grid)$dC +
                   inflow.grid-outflow.rate.grid*CONC

  # if(t-delay>0){
  #   #tmp<-lagvalue(t-delay)
  #   browser()
  # }

  return (list(dCONC))
  })
}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2, N2] <- ini  # initial concentration in the central point...
y.ini<- y

# solve for 8 time units
times <- 0:12
outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = parms,
               dim = c(N, N), lrw = 160000,CONCini=y.ini,delay=1,method= "lsodes")

outtimes <- seq(3,12,by=3)
image(outc, ask = FALSE, mfrow = c(2, 2),
      main = paste("time", outtimes),
      legend = TRUE, add.contour = TRUE,
      subset = time %in% outtimes)




Di erential Equations in R - R package deSolve<http://desolve.r-forge.r-project.org/slides/tutorial.pdf>
desolve.r-forge.r-project.org
Di erential Equations in R Tutorial useR conference 2011 Karline Soetaert, & Thomas Petzoldt Centre for Estuarine and Marine Ecology (CEME) Netherlands Institute of Ecology (NIOO-KNAW)



	[[alternative HTML version deleted]]