model <- function(t, state, parms) { with(as.list(c(state,parms)), { dR <- r*R*(1 - R/K) - a*(R^n)*N/(h^n+R^n) dN <- c*a*(R^n)*N/(h^n+R^n) - delta*N return(list(c(dR, dN))) }) } p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4,n=2) s <- c(R=1,N=0.01) plane() # Make Figures # size<-5#inch # # pdf("sigmoidA.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.485,n=2) # plane(xmax=1.5) # run(500,traject=TRUE) # newton(c(R=0.5,N=0.6),plot=TRUE) # newton(c(R=1,N=0),plot=TRUE) # newton(c(R=0,N=0),plot=TRUE) # dev.off() # # pdf("sigmoidB.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.42,n=2) # plane(legend=FALSE) # run(200,0.5,traject=TRUE) # newton(c(R=0.4,N=0.6),plot=TRUE) # newton(c(R=1,N=0),plot=TRUE) # newton(c(R=0,N=0),plot=TRUE) # dev.off() # # pdf("sigmoidC.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.15,n=2) # plane(legend=FALSE) # run(200,0.5,traject=TRUE) # newton(c(R=0.4,N=0.6),plot=TRUE) # newton(c(R=1,N=0),plot=TRUE) # newton(c(R=0,N=0),plot=TRUE) # dev.off() # # pdf("sigmoidD.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.4,a=1,c=1,delta=0.5,n=2) # plane(legend=FALSE) # run(100,0.5,traject=TRUE) # newton(c(R=0.4,N=0.6),plot=TRUE) # newton(c(R=1,N=0),plot=TRUE) # newton(c(R=0,N=0),plot=TRUE) # dev.off() # # pdf("sigmoidE.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.49,n=2) # f <- newton(c(R=0.7,N=0.4)) # continue(f,x="K",y="R",xmin=0.1,xmax=3,ymin=-0.1,ymax=2,step=0.005,positive=TRUE) # continue(c(R=1,N=0),x="K",y="R",xmin=0.1,xmax=3,ymin=-0.1,ymax=2,positive=TRUE,add=TRUE) # continue(c(R=0,N=0),x="K",y="R",xmin=0.1,xmax=3,ymin=-0.1,ymax=2,positive=TRUE,add=TRUE) # p["K"] <- 1.43 # Start at Hopf bifurcation # f <- newton(c(R=0.7,N=0.7)) # burnin <- 5000; tmax <- 500 # for (k in seq(1.43,3,0.01)) { # p["K"] <- k # 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(k,nrow(datamin)),datamin$R,pch=".") # points(rep(k,nrow(datamax)),datamax$R,pch=".") # } # dev.off() # # pdf("sigmoidF.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.49,n=2) # f <- newton(c(R=0.7,N=0.4)) # continue(f,x="K",y="N",xmin=0.1,xmax=3,ymin=-0.1,ymax=2,step=0.005,positive=TRUE) # continue(c(R=1,N=0),x="K",y="N",xmin=0.1,xmax=3,ymin=-0.1,ymax=2,positive=TRUE,add=TRUE) # p["K"] <- 1.43 # Start at Hopf bifurcation # f <- newton(c(R=0.7,N=0.7)) # burnin <- 5000; tmax <- 500 # for (k in seq(1.43,3,0.01)) { # p["K"] <- k # 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$N[2:(nr-1)]data$N[1:(nr-2)] & data$N[2:(nr-1)]>data$N[3:nr])),] # points(rep(k,nrow(datamin)),datamin$N,pch=".") # points(rep(k,nrow(datamax)),datamax$N,pch=".") # } # dev.off() # # pdf("sigmoidG.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.49,n=2) # f <- newton(c(R=0.7,N=0.4)) # continue(f,x="delta",y="N",xmin=0.01,xmax=0.55,ymin=-0.1,step=0.005,positive=TRUE) # continue(c(R=1,N=0),x="delta",y="N",xmin=0.01,xmax=0.55,ymin=-0.1,positive=TRUE,add=TRUE) # p["delta"] <- 0.28 # Start at Hopf bifurcation # f <- newton(c(R=0.2,N=0.3)) # burnin <- 1000; tmax <- 250 # for (k in seq(0.28,0.4775,0.005)) { # p["delta"] <- k # 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$N[2:(nr-1)]data$N[1:(nr-2)] & data$N[2:(nr-1)]>data$N[3:nr])),] # points(rep(k,nrow(datamin)),datamin$N,pch=".") # points(rep(k,nrow(datamax)),datamax$N,pch=".") # } # dev.off() # # pdf("sigmoidH.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.49,n=2) # p["K"] <-0.5; plane(xmax=3,ymax=1.8) # p["K"] <- 1; plane(xmax=3,ymax=1.8,add=TRUE,legend=FALSE);newton(c(R=0.75,N=1),plot=TRUE) # p["K"] <-1.5; plane(xmax=3,ymax=1.8,add=TRUE,legend=FALSE);newton(c(R=0.75,N=1),plot=TRUE) # p["K"] <- 2; plane(xmax=3,ymax=1.8,add=TRUE,legend=FALSE);newton(c(R=0.75,N=1),plot=TRUE) # p["K"] <-2.5; plane(xmax=3,ymax=1.8,add=TRUE,legend=FALSE);newton(c(R=0.75,N=1),plot=TRUE) # p["K"] <- 3; plane(xmax=3,ymax=1.8,add=TRUE,legend=FALSE); newton(c(R=0.75,N=1),plot=TRUE) # dev.off() # # pdf("sigmoidI.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=3,h=0.1,a=0.5,c=1,delta=0.49,n=2) # plane(xmax=3,ymax=1.8) # newton(c(R=0.75,N=1),plot=TRUE) # s["R"] <- 3 # run(1000,0.1,traject=TRUE) # dev.off() # # # pdf("sigmoidJ.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4,n=10) # s <- c(R=1,N=0.01) # plane(xmin=0,ymin=0,x="N",y="R",xmax=0.65,show="R",legend=FALSE) # dev.off() # # pdf("sigmoidK.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.485,n=2) # layout(matrix(c(1,2),ncol=1,nrow=2),heights=c(1,0.5)) # par(mar=c(0,0.1,0.1,0.1)) # plane() # par(mar=c(0.1,0.1,0,0.1)) # with(as.list(p),{curve(1-(R^n)/(h^n+R^n),from=0,to=1,xname="R",col="darkgreen",lwd=2)}) # legend("bottomleft",legend="f(R)",col="darkgreen",lty=1,lwd=2,cex=sizeLegend) # dev.off() # # pdf("sigmoidL.pdf",width=size,height=size) # par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i") # p <- c(r=1,K=1,h=0.4,a=1,c=1,delta=0.5,n=2) # layout(matrix(c(1,2),ncol=1,nrow=2),heights=c(1,0.5)) # par(mar=c(0,0.1,0.1,0.1)) # plane() # par(mar=c(0.1,0.1,0,0.1)) # with(as.list(p),{curve(1-(R^n)/(h^n+R^n),from=0,to=1,xname="R",col="darkgreen",lwd=2)}) # legend("bottomleft",legend="f(R)",col="darkgreen",lty=1,lwd=2,cex=sizeLegend) # dev.off()