################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap05") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=75, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) ## palette(gray((0:8)/8)) library(lattice) trellis.par.set(canonical.theme(color=FALSE)) myQ <- function(w=4,h=4, mr = c(5,4,.5,1.5)) { quartz(width=w,height=h) par(mar=mr) } ################################################### ### chunk number 2: weeds5 ################################################### #weeds <- read.csv("OldfieldPerennials.csv") library(primer) data(weeds) #trellis.device(color=F, width=5, height=5) print( xyplot(Cover~ Age, groups=Species, data=weeds, main="Oldfield Succession", #subset=Age < 20, lwd=2, ylim=c(2,13), xlim=c(2,21), #lty=c(1,1,3,2,4,5), lwd=c(1,2,1,2,1,1), xlab="Age (years)", scales=list(tck=c(1,0), y=list(at=seq(0,14,by=2) ), x=list(at=seq(0,25, by=5)) ), auto.key=list(points=F, columns=2, lines=T), # legend=list("inside"), type=c('smooth'), ylab="Percent Cover") ) ################################################### ### chunk number 3: ################################################### dlvcomp2 <- function(N, alpha, rd=c(1,1)) { N1.t1 <- N[1] + rd[1] * N[1] * (1 - alpha[1,1]*N[1] - alpha[1,2]*N[2]) N2.t1 <- N[2] + rd[2] * N[2] * (1 - alpha[2,1]*N[1] - alpha[2,2]*N[2]) c(N1.t1, N2.t1) } ################################################### ### chunk number 4: ################################################### alphs <- matrix(c( .01, .005, .008, .01), ncol=2, byrow=TRUE) t <- 20 ################################################### ### chunk number 5: ################################################### N <- matrix(NA, nrow=t+1, ncol=2) N[1,] <- c(10,10) for(i in 1:t) N[i+1,] <- dlvcomp2(N[i,], alphs) ################################################### ### chunk number 6: InterComp5 ################################################### matplot(0:t, N, type='l', col=1, ylim=c(0,110)) abline(h=1/alphs[1,1], lty=3) text(0, 1/alphs[1,1], "K", adj=c(0,0)) legend("right", c(expression("Sp.1 "*(alpha[21]==0.008)), expression("Sp.2 "*(alpha[12]==0.005))), lty=1:2, bty='n') ################################################### ### chunk number 7: ################################################### lvcomp2 <- function(t, n, parms) { with(as.list(parms), { dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2]) dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1]) list(c(dn1dt, dn2dt)) } ) } ################################################### ### chunk number 8: ################################################### library(deSolve) parms <- c(r1=1, r2=.1, a11=.2, a21=0.1, a22=0.02, a12=0.01); initialN <- c(2,1) out <- ode(y=initialN, times=1:100, func=lvcomp2, parms=parms) matplot(out[,1],out[,-1], type='l') ################################################### ### chunk number 9: ################################################### a <- matrix(c( 0.01, 0.005, 0.005, 0.01), ncol=2, byrow=TRUE) ################################################### ### chunk number 10: ################################################### N2iso <- expression(1/a[2,2] - (a[2,1]/a[2,2])*N1) ################################################### ### chunk number 11: ################################################### N1 <- 0:200 plot(N1, eval(N2iso), type='l', ylim=c(0, 200), xlim=c(0, 200), ylab=expression("N"[2])) ################################################### ### chunk number 12: ################################################### arrows(x0=90, y0=150, x1=90, y1=80, length=0.1) arrows(x0=75, y0=0, x1=75, y1=50, length=0.1) ################################################### ### chunk number 13: ################################################### a <- matrix(c( 0.01, 0.005, 0.005, 0.01), ncol=2, byrow=TRUE) myQ() curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, lty=2, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1]), pty="s" ) box() axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) arrows(c(75, 175), c(150,150), lty=2, c(75,175), c(80,45), length=0.12) arrows(c(25, 125), c(5, 5), lty=2, c(25,125), c(55, 20), length=0.12) dev.print(pdf, "N2iso.pdf") ### N1 isocline ### To graph this isocline, we have to solve the N1 isocline ### for N2 merely because N2 happens to be on the y axis myQ() #quartz(width=5.5/2.54, height=6.5/2.54) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1]),pty='s' ) axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) box() arrows(c(150,150), c(75, 175), c(80,45), c(75,175), length=0.12) arrows( c(5, 5), c(25, 125), c(55, 20), c(25,125), length=0.12) dev.print(pdf, "N1iso.pdf") ################################################### ### chunk number 14: ################################################### ### Isoclines for combinations of invasion criteria ### Both isoclines on the same graph a <- matrix(c( 0.01, 0.005, 0.005, 0.01), ncol=2, byrow=TRUE) myQ() #quartz(width=55/2.54, height=6.5/2.54) plot( (a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), (a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2, pch=1, col=1, axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s') curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T) axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) box() arrows(c(25,25,25), c(25,25,25), c(25,50,50), c(50,25,50), col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1)) arrows(c(150,150,150), c(150,150,150), c(150,120, 120), c(120,150,120), col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1)) arrows(c(10, 10, 10), c(125,125,125), c(10,30,30), c(105,125,105), length=0.12, col=1, lwd=c(1,1,2), lty=c(2,1,1)) arrows(c(125,125,125), c(10, 10,10), c(125,105,105), c(30,10,30), length=0.12, col=1, lwd=c(1,1,2), lty=c(2,1,1)) dev.print(dev=pdf, file="LVisocline1.pdf") ### Stable coexistence SEE ABOVE ### Species 2 can invade , species one cannot. myQ() a <- matrix(c( 0.01, 0.02, 0.005, 0.01), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2, axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T) axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), lwd=1, c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) box() xs <- c(1, 1, .4, .4) * a[1,1]^-1 ys <- c(.2, .4, .4, .7)*a[2,2]^-1 arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4], lty=2:1, lwd=c(1,1), length=0.12) points(0,1/a[2,2], cex=2) dev.print(dev=pdf, file="LVisocline2wins.pdf") ### Species 1 can invade , species 2 cannot. myQ() a <- matrix(c( 0.01, 0.005, 0.02, 0.01), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2, axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T) box() axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) ys <- c(1, 1, .4, .4) * a[1,1]^-1 xs <- c(.2, .4, .4, .7)*a[2,2]^-1 arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4], lty=1:2, lwd=c(1,1), length=0.12) points(1/a[1,1],0, cex=2) dev.print(dev=pdf, file="LVisocline1wins.pdf") ### Neither species can invade , yet an unstable equilbrium -- a saddle point attractor-repellor -- exists. myQ() a <- matrix(c( 0.01, 0.02, 0.02, 0.01), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lty=2, axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T) box() axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) pxys <- matrix(c(.13,.1,.25,.1, .1,.13,.1,.25, .115,.115,.25,.25, .7,.8,.4,.8, .8,.7,.8,.4, .75,.75,.4,.4, .25,.4,.1, .6, .4,.25,.6, .1), nc=4, byrow=T) xys <- pxys xys[,1:2] <- a[1,1]^-1 * pxys[,1:2] xys[,3:4] <- pxys[,3:4] * a[2,2]^-1 arrows(xys[,1], xys[,2], xys[,3], xys[,4], lty=c(1,2,3,1,2,3,2,1), lwd=c(1), length=0.12) n2star <- (a[1,1] - a[2,1])/(a[1,1]*a[2,2] - a[1,2]*a[2,1]) n1star <- (a[2,2] - a[1,2])/(a[1,1]*a[2,2] - a[1,2]*a[2,1]) points(n1star,n2star, cex=2) dev.print(dev=pdf, file="LVisocline0win.pdf") ################################################### ### chunk number 15: ################################################### N1Star <- expression( (a22-a12)/(a22*a11 - a12*a21) ) N2Star <- expression( (a11-a21)/(a22*a11 - a12*a21) ) ################################################### ### chunk number 16: ################################################### a11 <- a22 <- 0.01; a12 <- 0.001; a21 <- 0.001 N1 <- eval(N1Star) N2 <- eval(N2Star) N1 ################################################### ### chunk number 17: ################################################### dN1dt <- expression( r1 * N1 - r1 * a11 * N1^2 - r1 * a12 * N1 * N2 ) dN2dt <- expression(r2*N2 - r2*a22*N2^2 - r2*a21*N1*N2) ################################################### ### chunk number 18: ################################################### ddN1dN1 <-D(dN1dt, "N1") ddN1dN1 ################################################### ### chunk number 19: ################################################### ddN1dN2 <- D(dN1dt, "N2") ddN2dN1 <- D(dN2dt, "N1") ddN2dN2 <- D(dN2dt, "N2") ################################################### ### chunk number 20: ################################################### J <- expression( matrix(c(eval(ddN1dN1), eval(ddN1dN2), eval(ddN2dN1), eval(ddN2dN2)), nrow=2, byrow=TRUE) ) ################################################### ### chunk number 21: ################################################### r1 <- r2 <- 1 J1 <- eval(J) J1 ################################################### ### chunk number 22: ################################################### eigStable <- eigen(J1); eigStable[["values"]] ################################################### ### chunk number 23: ################################################### a11 <- a22 <- 0.01 a12 <- a21 <- 0.011 N1 <- eval(N1Star); N2 <- eval(N2Star) eigen( eval(J) )$values ################################################### ### chunk number 24: ################################################### a11 <- a21 <- 0.01; a22 <- a12 <- 0.015 ################################################### ### chunk number 25: ################################################### N1 <- N2 <- 1/(a11+a22) eigen( eval(J) )[["values"]] ################################################### ### chunk number 26: ################################################### N1 <- 1/(a11); N2<-0 eigen( eval(J) )[["values"]] ################################################### ### chunk number 27: ################################################### ### Rates of change with respect to each other: partials N1Star <- expression( (a22-a12)/(a22*a11 - a12*a21) ) N2Star <- expression( (a11-a21)/(a22*a11 - a12*a21) ) dN1dt <- expression(r1*N1 - r1*a11*N1^2 - r1*a12*N1*N2) dN2dt <- expression(r2*N2 - r2*a22*N2^2 - r2*a21*N1*N2) ### respect to each other. (ddN1dN1 <- D(dN1dt, "N1")) (ddN1dN2 <- D(dN1dt, "N2")) (ddN2dN1 <- D(dN2dt, "N1")) (ddN2dN2 <- D(dN2dt, "N2")) ### Tell R what the alpha's are a <- matrix(c( 0.01, 0.00, 0.005, 0.01), ncol=2, byrow=TRUE) r1 <- r2 <- r <- 1 a11 <- a[1,1]; a12 <- a[1,2]; a21 <- a[2,1]; a22 <- a[2,2] ### Evaluate the Jacobian using the above parameters N1 <- eval(N1Star); N2 <- eval(N2Star) ddN1dN1 <- D(dN1dt, "N1") ddN1dN2 <- D(dN1dt, "N2") ddN2dN1 <- D(dN2dt, "N1") ddN2dN2 <- D(dN2dt, "N2") J <- matrix(c(eval(ddN1dN1), eval(ddN1dN2), eval(ddN2dN1), eval(ddN2dN2)), nrow=2, byrow=TRUE) J (ea <- eigen(J)) LSA <- function(N1Star.exp=N1Star, N2Star.exp=N2Star, dN1dt.exp=dN1dt, dN2dt.exp=dN2dt, r1=r, r2=r, aM=a) { a11 <- aM[1,1]; a12 <- aM[1,2]; a21 <- aM[2,1]; a22 <- aM[2,2] ### Evaluate the Jacobian using the above parameters N1 <- eval(N1Star.exp); N2 <- eval(N2Star.exp) ddN1dN1 <- D(dN1dt.exp, "N1") ddN1dN2 <- D(dN1dt.exp, "N2") ddN2dN1 <- D(dN2dt.exp, "N1") ddN2dN2 <- D(dN2dt.exp, "N2") J <- matrix(c(eval(ddN1dN1), eval(ddN1dN2), eval(ddN2dN1), eval(ddN2dN2)), nrow=2, byrow=TRUE) ea <- eigen(J) list(r.s=c(r1,r2), a=aM, Jacobian=J, e.vals =ea$values, e.vecs<-ea$vectors, pdes = list(ddN1dN1,ddN1dN2, ddN2dN1, ddN2dN2) ) } library(deSolve) #source("LVComp2.r") params <- c(r1=1, r2=1, a11 = a[1,1], a12 = a[1,2], a21 = a[2,1], a22 = a[2,2]) a12s <- seq(.0041, .0161, length=100) lambdas <- t( sapply(as.list(a12s), function(x) {a[1,2] <- x; a[2,1] <- x LSA(aM=a)$e.vals}) ) myQ() matplot(a12s/params["a11"], lambdas, type="l", ylab="Eigenvalues", xlab=expression(italic( alpha["ij"]/alpha["ii"] )), lty=1 ) abline(h=0, lty=3) abline(v=1, lty=3) text(.6, -.33, "Stable", adj=0 ) text(1.4, .05, "Unstable" ) dev.print(dev=pdf, file="LVStab.pdf") myQ() a <- matrix(c( 0.01, 0.005, 0.005, 0.01), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T, lwd=2, lty=2) axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) box() nstart <- 10 n1 <- c(seq(1/(2*a[1,1]), 1/a[2,1], length=nstart), seq(3/(4*a[1,1]), 1, length=nstart) ) n2 <- c(seq(1/a[1,2], 1/(2*a[2,2]), length=nstart), seq(1, 3/(4*a[2,2]), length=nstart ) ) paramsST <- c(r1=.1, r2=.1, a11 = a[1,1], a12 = a[1,2], a21 = a[2,1], a22 = a[2,2]) for( i in 1:length(n1) ) {Ns <- ode(c(n1[i],n2[i]), 1:200, lvcomp2, paramsST ) points(n1[i],n2[i], pch=19) lines(Ns[,2],Ns[,3]) } dev.print(dev=pdf, file="LVStable.pdf") ### Dynamics of Saddle repellors ### Neither species can invade , yet an unstable equilbrium -- a saddle point attractor-repellor -- exists. myQ() a <- matrix(c( 0.01, 0.02, 0.02, 0.01), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T, lwd=2,lty=2) axis(1, c(0, 1/a[1,1], 1/a[2,1]), c(0, expression(1/alpha[11]),expression(1/alpha[21])) ) axis(2, c(0, 1/a[2,2], 1/a[1,2]), c(0, expression(1/alpha[22]),expression(1/alpha[12])) ) box() nstart <- 10 n1 <- c(seq(1/(2*a[1,1]), 1/(a[1,1]), length=nstart), seq(1/(2*a[2,1]), 0, length=nstart) ) n2 <- c(seq(1/a[2,2], 1/(2*a[2,2]), length=nstart), seq(0, 1/(2*a[1,2]), length=nstart ) ) paramsUN <- c(r1=.1, r2=.1, a11 = a[1,1], a12 = a[1,2], a21 = a[2,1], a22 = a[2,2]) for( i in 1:length(n1) ) {Ns <- ode(c(n1[i],n2[i]), 1:100, lvcomp2, paramsUN ) points(n1[i],n2[i], pch=19) lines(Ns[,2],Ns[,3])} #abline(c(0,1), lwd=2, col="forestgreen") dev.print(dev=pdf, file="LVSaddle.pdf") myQ() a <- matrix(c( 0.01, 0.015, 0.01, 0.015), ncol=2, byrow=TRUE) curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], axes=FALSE, ylim=c(0,100), xlim=c(0, 100), ylab=expression("N"[2]), xlab=expression("N"[1])) curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T, lwd=4, lty=2) axis(1, c(0, 1/a[1,1]), c(0, expression(1/alpha[11])) ) axis(2, c(0, 1/a[2,2]), c(0, expression(1/alpha[22])) ) box() nstart <- 6 n1 <- c(seq(0, 1/(a[1,1]), length=nstart), seq(1/(2*a[1,1]), 0, length=nstart) ) n2 <- c(seq(1.5/a[2,2], 1/(2*a[2,2]), length=nstart), seq(0, 1/(2*a[2,2]), length=nstart ) ) paramsUN <- c(r1=1, r2=1, a11 = a[1,1], a12 = a[1,2], a21 = a[2,1], a22 = a[2,2]) for( i in 2:length(n1) ) {Ns <- ode(c(n1[i],n2[i]), 1:10, lvcomp2, paramsUN ) points(n1[i],n2[i], pch=19) lines(Ns[,2],Ns[,3]) arrows(Ns[1,2],Ns[1,3],Ns[10,2],Ns[10,3], length=0.12)} points(1/sum(a[1,]), 1/sum(a[1,]), cex=2) dev.print(dev=pdf, file="LVNeutral.pdf") ratio <- ode(c(1,10), 1:100, lvcomp2, paramsUN );ratio[,2]/ratio[,3] a ################################################### ### chunk number 28: ################################################### findL1 <- function(B, r1=1, r2=1, a11=.01, a22=.01) { r1=r1; r2=r2; a11 = a11; a22 <- a22 a12 <- B[1]*a22; a21 <- B[2]*a11 if(B[1] < 1 & B[2] > 1) {N1 <- 1/a11 N2=0 } else if(B[1] > 1 & B[2] < 1) {N1 <- 0; N2 <- 1/a22 } else if(B[1] < 1 & B[2] < 1) { N1 <- eval(N1Star); N2=eval(N2Star) } else if(B[1] > 1 & B[2] > 1) {N1 <- eval(N1Star) N2 <- eval(N2Star) } else ### This is completely arbitrary if(B[1] == 1 & B[2] == 1) N1 <- N2 <- 1/(2*a11) ### Evaluate the Jacobian using the above parameters J <- expression( matrix(c(eval(ddN1dN1), eval(ddN1dN2), eval(ddN2dN1), eval(ddN2dN2)), nrow=2, byrow=TRUE) ) J1 <- eval(J); J1 eigAll <- eigen(J1) return( list(domL=max(Re(eigAll$values)), eigenanal = eigAll, Jacobian=J1, N = c(N1, N2) ) ) } a12/a22 findL1(B=c(1.5,1/1.5)) ### Figs for Stable Eq dat <- expand.grid(B12=seq(0.0001, 1, by=0.05), B21=seq(0.0001, 1, by=0.05)) dat$L <- apply(dat, 1, function(b) findL1(B=b)$domL) trellis.device(col=FALSE, width=4, height=4) print(wireframe(L ~ B12 + B21, dat, scales=list(arrows=FALSE), screen = list(z = 20, x = -70), xlim=c(0,1), ylim=c(0,1), drape=TRUE, colorkey=TRUE, zlab=expression(lambda[1]), xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeig3D.pdf") trellis.device(col=FALSE, width=4, height=4) print(contourplot(L ~ B12 + B21, dat, cuts=10, xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeigContour.pdf") ### Figs for unstable eq dat2 <- expand.grid(B12=seq(1.01, 2.01, by=0.05), B21=seq(1.01, 2.01, by=0.05)) dat2$L <- apply(dat2, 1, function(b) findL1(B=b)$domL) trellis.device(col=FALSE, width=4, height=4) print(wireframe(L ~ B12 + B21, dat2, scales=list(arrows=FALSE), screen = list(z = 20, x = -70), xlim=c(1,2), ylim=c(1,2), drape=TRUE, colorkey=T, zlab=expression(lambda[1]), xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeig3Du.pdf") trellis.device(col=FALSE, width=4, height=4) print(contourplot(L ~ B12 + B21, dat2, cuts=10, ylim=c(.95,2.05), xlim=c(.95, 2.05), xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeigContouru.pdf") ### Games dat3 <- expand.grid(B12=seq(.01, 2, by=0.05), B21=seq(.01, 2, by=0.05)) dat3$one <- dat3$B12*dat3$B21 dat3$L <- apply(dat3, 1, function(b) findL1(B=b)$domL) dat3$RT <- dat3$L^(-1) trellis.device(col=FALSE, width=4, height=4) print(contourplot(L ~ B12 + B21, dat3, cuts=20, labels=T, xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeigContourAll.pdf") print(wireframe(L ~ B12 + B21, dat3, scales=list(arrows=FALSE), screen = list(z = 20, x = -70), zlim=c(-1,.35), drape=TRUE, colorkey=F, zlab=expression(lambda[1]), xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeig3DAll.pdf") dat3$L2 <- apply(dat3, 1, function(b) findL1(B=b, r1=.5, r2=.5)$domL) trellis.device(col=FALSE, width=4, height=4) print(contourplot(L2 ~ B12 + B21, dat3, cuts=20, labels=T, xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeigContourAllhalf.pdf") print(wireframe(L2 ~ B12 + B21, dat3, scales=list(arrows=FALSE), screen = list(z = 20, x = -70), zlim=c(-1,.35), drape=TRUE, colorkey=F, zlab=expression(lambda[1]), xlab=expression(beta[12]), ylab=expression(beta[21]))) dev.print(pdf, "LVeig3DAllhalf.pdf")