# Filename: LeslieFecundity.R # R script to # - Find the dominant eigenvalue and corresponding # eigenvectors for an array of juvenile fecundities # - Print out a table of the results # - Print out results graphically # Print out the header for the table cat("Juvenile Dominant Proportion Structure \n") cat("Fecundity Eigenvalue Hatchlings Juveniles Adults\n") cat("------------------------------------------------------\n") # Array of juvenile fecundity values f = c(0.1, 0.2, 0.5, 1.5, 2, 2.5, 3, 4) # Make a vector that we will use later normv = t(t(numeric(3))) # Start the loop for (i in 1:length(f)){ # loop through length of f vector # Construct Leslie matrix with appropriate juvenile fecundity A = matrix( c( 0, 0.5, 0, f[i], 0, 0.25, 3, 0, 0 ), ncol = 3) # Get eigenvalues and eigenvectors lambda = Re(eigen(A)$values) v = Re(eigen(A)$vectors) # Find dominant eigenvalue in lambda # and make note of position in lambda matrix lambdavector = abs(lambda) # array of eigenvalues for (j in 1:length(lambdavector)){ if (max(abs(lambda)) == abs(lambdavector[j])){ loc = j # position of dom eig in lambda vector deig = lambda[j] # domianant eigenvalue } } # Normalize eigenvector associated with dominant eigenvalue normv = cbind(normv, v[,loc]/sum(v[,loc])) # This creates a matrix called normv in which the # i^th column corresponds to the dominant eigenvector # associated with the i^th juvenile fecundity value # Print these values cat(sprintf(" %3.1f ",f[i])) # juvenile fecundity value cat(sprintf(" %5.3f ",deig)) # dominant eigenvalue cat(sprintf(" %5.3f ",normv[1,i+1])) # hatchling prop. at eq cat(sprintf(" %5.3f ",normv[2,i+1])) # juv. prop. at eq cat(sprintf(" %5.3f ",normv[3,i+1])) # adult prop. at eq cat("\n") # new line } # Remove the first column from normv normv = normv[,-1] # Generate graph h = normv[1,] # vector of hatchling proportions at eq j = normv[2,] # vector of juvenile proportions at eq a = normv[3,] # vector of adult proportions at eq plot(f, h, col = "red", type = "b", pch = 20, xlab = "Juvenile Fecundity", ylab = "Equilibrium Structure", xlim = c(0, 4), ylim = c(0, 1)) par(new = TRUE) plot(f, j, col = "green", type = "b", pch = 20, xlab = "", ylab = "", xaxt = "n", yaxt = "n", xlim = c(0, 4), ylim = c(0, 1)) par(new = TRUE) plot(f, a, col = "blue", type = "b", pch = 20, xlab = "", ylab = "", xaxt = "n", yaxt = "n", xlim = c(0, 4), ylim = c(0, 1)) legend("topright", c("Hatchlings", "Juveniles", "Adults"), col = c("red", "green", "blue"), lty = c(1,1,1), pch = c(20, 20, 20))