numerical example of next generation matrix

NGS.R

iter_nextstep <- function(n, K){
  I <- c(1, 0)
  I_hist <- matrix(0, nrow = 2, ncol = n+1)
  I_hist[, 1] <- I
  for(i in 1:n){
    I <- K%*%I
    I_hist[,i+1] <- I  
  }
  return(t(I_hist))  
}

K <- matrix(c(1.4, 0.4,
              0.4, 0.3), byrow = TRUE, nrow = 2)
ei_K <- eigen(K)
norm_v <- ei_K$vectors[,1]/sum(ei_K$vectors[,1])

n <- 10
I_hist <- iter_nextstep(n,K)

matplot(0:n, I_hist, type = "l")
sumI <- rowSums(I_hist)
normI <- sweep(I_hist,1,sumI,"/")

png("fig1.png")
plot(sumI[-1]/sumI[-length(sumI)],
     ylab = expression(N[t+1]/N[t]), xlab = "t")
abline(h=ei_K$values[1],lty=3)
dev.off()

png("fig2.png")
matplot(0:n, normI, type = "p",ylab = expression(I[t]/N[t]), xlab = "t")
abline(h = norm_v, lty=3)
dev.off()