model <- function(t, state, parms) { with(as.list(c(state,parms)), { dR <- r*R*(1 - R/K) - a*R*N/(h+R) dN <- c*a*R*N/(h+R) - delta*N return(list(c(dR, dN))) }) } p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4) s <- c(R=1,N=0.01) plane() newton(c(R=0,N=0),plot=TRUE) newton(c(R=1,N=0),plot=TRUE) newton(c(R=0.4,N=0.6),plot=TRUE) # Make Figures # size<-5#inch # # pdf("monodA.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.46) # plane(xmax=1.5) # newton(c(R=1,N=0),plot=TRUE) # newton(c(R=0,N=0),plot=TRUE) # dev.off() # # pdf("monodB.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) # plane(legend=FALSE) # run(200,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("monodC.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) # plane(legend=FALSE) # run(200,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("monodD.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=1,a=1,c=1,delta=0.4) # plane(legend=FALSE) # run(200,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("monodE.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) # f <- newton(c(R=0.4,N=0.6)) # continue(f,x="K",y="R",xmin=0.1,xmax=1.5,ymin=-0.1,step=0.005,positive=TRUE) # continue(c(R=1,N=0),x="K",y="R",xmin=0.1,xmax=1.5,ymin=-0.1,positive=TRUE,add=TRUE) # continue(c(R=0,N=0),x="K",y="R",xmin=0.1,xmax=1.5,ymin=-0.1,positive=TRUE,add=TRUE) # p["K"] <- 0.9 # Start at Hopf bifurcation # f <- newton(c(R=0.4,N=0.6)) # burnin <- 1e4; tmax <- 500 # for (k in seq(0.9,1.5,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("monodF.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) # f <- newton(c(R=0.4,N=0.6)) # continue(f,x="K",y="N",xmin=0.1,xmax=1.5,ymin=-0.1,step=0.005,positive=TRUE) # continue(c(R=1,N=0),x="K",y="N",xmin=0.1,xmax=1.5,ymin=-0.1,positive=TRUE,add=TRUE) # p["K"] <- 0.9 # Start at Hopf bifurcation # f <- newton(c(R=0.4,N=0.6)) # burnin <- 1e4; tmax <- 500 # for (k in seq(0.9,1.5,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("monodG.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) # run(300,legend=FALSE) # dev.off()