Theta-logistic comparison

The following figure:

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()
```

page revision: 3, last edited: 13 May 2008 22:03