[R-sig-dyn-mod] better implementation of variable species ode.2D model

[Thomas Petzoldt]

Hi,

loops are not always bad, but one should avoid them in functions that 
are repeatedly called. As said, matrix-computations can help. See the 2D 
L&V model at page 10 of

http://dx.doi.org/10.18637/jss.v033.i09

or the 2D SIR model below with or without ReacTran.

ThPe



-------------- next part --------------
library(deSolve)
library(ReacTran)


## The model equations
SIR2D <- function (t, y, parms)  {
  S  <- matrix(y[ndx], N, N)
  I  <- matrix(y[N^2 + ndx], N, N)
  R  <- matrix(y[2*(N^2) + ndx], N, N)

  infect <- beta * I * S
  recovr <- nu * I

  dS <- -infect
  dI <- infect - recovr
  dR <- recovr

  if (reactran) {
    dI <- dI + tran.2D(I, dx = dx, dy = dy, D.x = D, D.y = D)$dC
  } else {
    Flux.Ix <- -D * rbind(zero, (I[2:N,] - I[1:(N-1),]), zero)/dx
    Flux.Iy <- -D * cbind(zero, (I[,2:N] - I[,1:(N-1)]), zero)/dy

    dFlux.I  <- (Flux.Ix[1:N,] - Flux.Ix[2:(N+1),])/dx +
                (Flux.Iy[,1:N] - Flux.Iy[,2:(N+1)])/dy

    dI <- dI + dFlux.I
  }
  list(c(dS, dI, dR))
}

## grid
N     <- 101
zero <- rep(0, N)
ndx <- 1:N^2
dx <- 10/N
dy <- 10/N

## model parameters
D     <- 1e-3
S0    <- 100
I0    <- 1
R0    <- 0.0
beta  <- 0.001
nu    <- 0.02


# initial conditions
S  <- matrix(nrow = N, ncol = N, data = S0)
I  <- matrix(nrow = N, ncol = N, data = 0)
R  <- matrix(nrow = N, ncol = N, data = R0)

set.seed(123)
I[sample(1:length(I), 10)] <- I0
#dim(I) <- c(N, N)  # redundant

y <- c(S = S, I = I, R = R)

times <- 0:200

reactran <- TRUE

system.time(
  out <- ode.2D (y = y, func = SIR2D, t = times, nspec=3, parms = NULL,
                dim = c(N, N), method = "adams")
)


windows(width=15, height=5)
image(out, select=1:3, legend = TRUE, ask=FALSE, main=c("Susceptible", "Infectious", "Recovered"), mfrow=c(1,3))