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

[Thomas Petzoldt]

Hi,

you are right that lsodes() supports DDEs, but ode.2D() does not and 
this has good reasons.

I don't know if it will work for you, but you may try to implement your 
model directly with dede() This requires of course that you set the 
appropriate structure of the Jacobian yourself. Note also that the "s" 
in lsodes stands for a "sparse" Jacobian and not for "stiff" equations. 
So you may also try lsoda etc.

BTW: I am not a great fan od DDEs, so I would prefer to employ 
intermediate pools to implement a time delayed system.

Thomas

On 01.11.2018 at 01:04 Raffaello Vardavas wrote:
> 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)