[R-sig-dyn-mod] Adding a delay to a ode.2d model
Thomas Petzoldt
thom@@@petzoldt @ending from tu-dre@den@de
Thu Nov 1 23:10:44 CET 2018
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)
More information about the R-sig-dynamic-models
mailing list