Skip to content
Prev 468 / 696 Next

[R-sig-dyn-mod] R-sig-dynamic-models Digest, Vol 106, Issue 1

Good day Jose,
Here is the code:

rm(list=ls())
library(ReacTran)
N <-200
alpha12<-0.0034
alpha21<-0.05
omega<-0.7
sigma<-1.2
rho<-1.5
mu<-0.02
gamma<-0.01
Grid <- setup.grid.1D(x.up = 0, x.down = 1, N = N)
Xini<-rep(x= 1, times= N)
Yini <- rep(x = 3, times = N)
Zini<-rep(x=4, times=N)
yini <- c(Xini, Yini, Zini)
Predator_prey <- function(t, y, parms) {
  X <- y[1:N]
  Y <- y[(N+1):(2*N)]
  Z <- y[(2*N+1):(3*N)]
  dX <- X*(1-X-alpha12*Y)-omega*X^2*Z/(1+X^2)+
    tran.1D (C = X, C.up = 0, C.down = 1,
             D = 0.1, dx = Grid)$dC+tran.1D (C = Z, C.up = NULL, C.down =
NULL,
                                              D = 1, dx = Grid)$dC
  dY <-  rho*Y*(1-Y-alpha21*X)+
    tran.1D (C = Y, C.up = NULL, C.down = NULL,
             D = 0.1, dx = Grid)$dC
  dZ <- sigma*X^2*Z/(1+X^2)-mu*Z-gamma*Z
    tran.1D (C = Z, C.up = NULL, C.down = NULL, flux.up = NULL, flux.down =
NULL,
             D = 0.1, dx = Grid)$dC+tran.1D (C = X, C.up = NULL, C.down =
NULL,
                                              D= -1, dx = Grid)$dC
  list(c(dX, dY, dZ))
}
times <- seq(from = 0, to = 100, by = 0.01)
print(system.time(
  out <- ode.1D(y = yini, func = Predator_prey,
                times = times, parms = NULL, nspec = 3,
                names = c("X", "Y","Z"), dimens = N)
))
matplot(out[,"time"], out[,2], type = "l", xlab = "time", ylab = "",
         main = "", lwd = 2)

Regards

Titus Kassem

On Thu, Jul 7, 2016 at 11:00 AM, <r-sig-dynamic-models-request at r-project.org

  
  
Message-ID: <CAG92UYAztQmJnOgyJSgUZRtnifiLtr-qX-+Zjwn9aUb-k=sirw@mail.gmail.com>
In-Reply-To: <mailman.17.1467885603.11670.r-sig-dynamic-models@r-project.org>