Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment