Instantly share code, notes, and snippets.

# tpoisot/isoclines.r Created Mar 18, 2014

Show isoclines and vector field
 library(RColorBrewer) library(shape) dndt = function(n1, n2, params) { with(as.list(params),{ dn1 = r1 * (1 - (a11 * n1 + a12 * n2)/k1) dn2 = r2 * (1 - (a22 * n2 + a21 * n1)/k2) return(list(dn1 = dn1, dn2 = dn2)) }) } params = list(r1 = 0.5, r2 = 0.5, a11 = 1.0, a12 = 0.6, a21 = 0.6, a22 = 1.0, k1 = 8, k2 = 8) step = 0.1 start = 0 end = 10 n1 = seq(from=start+step/2, to=end-step/2, by=step) n2 = seq(from=start+step/2, to=end-step/2, by=step) S = matrix(0, ncol=length(n1), nrow=length(n2)) for(i in c(1:length(n1))) { for(j in c(1:length(n2))) { dn = dndt(n1[i], n2[j], params) S[i,j] = sqrt(sum(unlist(dn))^2) } } im_col = colorRampPalette(brewer.pal(11, 'YlOrRd'))(150) par(pty='s') image(n1, n2, S, col=im_col, xlab='Species 1', ylab='Species 2') isoclines = function(params) { with(as.list(params),{ segments(0,k1/a12,k1,0, lwd=2, lty=1) segments(0,k2,k2/a21,0, lwd=2, lty=2) legend('topright', lty=c(1,2), lwd=2, legend=c('Sp. 1', 'Sp. 2')) }) } step = 0.5 start = 0 end = 10 n1 = seq(from=start+step/2, to=end-step/2, by=step) n2 = seq(from=start+step/2, to=end-step/2, by=step) for(i in c(1:length(n1))) { for(j in c(1:length(n2))) { dn = dndt(n1[i], n2[j], params) Arrows(n1[i], n2[j], n1[i]+dn\$dn1*0.01, n2[j]+dn\$dn2*0.01, arr.length=0.15, arr.type='triangle', segment=F) } } isoclines(params) box()