# Chaos in a three-species food chain: # Hastings & Powell (1991) Ecology 72. 896-903. # This illustrates: # 3-dimensional nullclines and chaotic trajectories # Takens reconstruction by shifting the numerical output # Effect of noise on a parameter or a variable on the attractor # Bifurcation diagram using continue() and by following # a limit cycle by tracking its minimum and maximum values. # Source both grind.R and cube.R! model <- function(t, state, parms) { with(as.list(c(state,parms)), { dR <- R*(1 - R/K) - c1*N*R/(1+b1*R) dN <- -aN*N + c1*N*R/(1+b1*R) - c2*Z*N/(1+b2*N) dZ <- -aT*Z + c2*Z*N/(1+b2*N) return(list(c(dR, dN, dZ))) }) } opar <- par() # Save graphics parameters par(mar=c(3.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)) p <- c(b1=6,b2=2,c1=5,c2=0.1,aN=0.4,aT=0.01,K=1) s <- c(R=1,N=0.1,Z=0.01) plane(eps=-.01) plane(eps=-.01,show=names(s),zero=F) f <- newton(c(R=0.1,N=0.23,Z=0),plot=T) run(50,tstep=0.1,traject=TRUE) s <- c(R=1,N=0.1,Z=0.1) f <- run(1e4) # Burn in cube(x=1,y=3,z=2,ymax=6,eps=1e-9) run3d(x=1,y=3,z=2,tmax=5000,tstep=0.2,state=f,add=TRUE) plotdev(theta=10,phi=20) # Plot same trajectory using the rgl library library(rgl) data <- run(tmax=5000,tstep=0.2,state=f,timeplot=FALSE,table=TRUE) plot3d(x=data$R,y=data$N,z=data$Z,xlab="R",ylab="N",zlab="Z",type="l") # Make a Takens reconstruction: data <- run(tmax=5000,state=f,table=TRUE) plot(data$R[1:4999],data$R[2:5000],pch=".") # Noise on a parameter, run and a Takens reconstruction: noiseP="parms[\"K\"]<-abs(rnorm(1,1,0.1))" data <- run(1e4,state=f,table=TRUE,after=noiseP) plot(data$R[5001:9999],data$R[5002:1e4],pch=".") # Noise on a variable, run and a Takens reconstruction: noiseS="if(runif(1)<0.1)state[3]=state[3]+abs(rnorm(1,0,0.1))" data <- run(1e4,state=f,table=TRUE,after=noiseS) plot(data$R[5001:9999],data$R[5002:1e4],pch=".") # Make a bifurcation diagram p["b1"] <- 1 s <- c(R=1,N=0.1,Z=0.1) f <- run(1e4) # Burn in continue(newton(f),x="b1",xmin=0,xmax=4,y="R",atol=1e-9,rtol=1e-9) p["b1"] <- 2.14 # Start at Hopf bifurcation burnin <- 1e4; tmax <- 500 for (b1 in seq(2.14,4,0.01)) { p["b1"] <- b1 f <- run(burnin,state=f,timeplot=FALSE) data <- run(tmax,tstep=0.2,state=f,table=TRUE,timeplot=FALSE) nr <- nrow(data) datamin <- data[(1+which(data$R[2:(nr-1)]data$R[1:(nr-2)] & data$R[2:(nr-1)]>data$R[3:nr])),] points(rep(b1,nrow(datamin)),datamin$R,pch=".") points(rep(b1,nrow(datamax)),datamax$R,pch=".") }