Theta-logistic comparison

The following figure:

thetalogist.png

Is produced by this code:

## closed-form solution for theta-logistic
tlogfun = function(t,theta=1.1,n0=0.1,r=1,K=1) {
    (K^(-theta)+(n0^(-theta)-K^(-theta))*exp(-r*theta*t))^(-1/theta)
}

## theta-logistic derivative function, using with() trick
##   to get variable and parameter names
tlderiv = function(time,y,parms) {
    L <- c(as.list(y),as.list(parms))
    gr <- with(L,r*n*(1-(n/K)^theta))
    list(gr,NULL)
}

## numeric solution using LSODA from the odesolve package
tlogderiv = function(t,theta=1.1,n0=0.1,r=1,K=1) {
    L = lsoda(y=c(n=n0),times=t,func=tlderiv,
      parms=c(r=r,K=K,theta=theta))
    L[,"n"]
}

library(odesolve) ## load package
thetavec = c(0.1,0.5,0.75,1,1.33,2,10)  ## define theta values
tvec = seq(0,10,by=0.2)                 ## times for solution
## calculate closed-form solutions for each theta
cform <- sapply(thetavec,
                tlogfun,t=tvec)
## calculate numeric solutions for each theta
lform <- sapply(thetavec,
                tlogderiv,t=tvec)

## doing clever stuff with colors -- sequences
##  from red to black and from black to blue
cols = c(hsv(h=1,s=1,v=seq(from=1,to=0,length=4)),
  hsv(h=0.5,s=1,v=seq(from=0,to=1,length=4)[-1]))

## set up to output to a PNG file (for web view)
png(width=1000,height=600,file="thetalogist.png")
## two graphs on one page, horizontal labels, squeeze margins
par(mfrow=c(1,2),las=1,mar=c(5,4,1,1)+0.1)
## base curve of dN/dt as a function of N (logistic)
curve((1-x)^1,from=0,to=1,xlab="pop density",ylab="per cap growth")
## add curves for different theta values
for (i in seq(along=thetavec)) {
    curve((1-x)^thetavec[i],add=TRUE,col=cols[i])
}
## plot theoretical curves
matplot(tvec,cform,type="l",lty=1,col=cols,
        xlab="pop density",ylab="time")
## add points for numerical solutions
matpoints(tvec,lform,pch=2,col=cols)
## add legend
legend(7,0.8,
       title=expression(theta),
       thetavec,col=cols,lty=1)
dev.off()
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License