Why the eigen function returns "number of items to replace is not a multiple of replacement length"?

I'm using for the question here the classic example of the Lotka-Volterra model with predation (2 ODEs, 6 parameters). I need to calculate the equilibrium point, which I know the analytical expression of, and the eigenvalues of the Jacobian matrix at this point. I need to do this for a big number of locations (here I'll do only 2), which differ in this example by the value of 2 parameters, epsilon and psi (2 parameter values for each for the 2 locations).

I created a loop (because again the size of epsilon and psi is going to be much bigger). Here is my code:

A21 = as.matrix(c(0, 0))
alpha = -1
beta = 2
gamma = 1
delta = -2
epsilon = as.matrix(c(0, 1))
psi = as.matrix(c(0, -2))
x = 0
y = 0
param = c(0,0,0,0,0,0)
eig = A21
eqn <- function (t, state, pars)
{
  with (as.list(c(state, pars)),  {
    dx <- x*(alpha - beta*y - epsilon)
    dy <- -y*(gamma - delta*x + psi)
    list(c(dx, dy))
  })
}

for(i in 1:dim(A21)[1]) {
  x[i] = (gamma + psi[i]) / delta
  y[i] = (alpha - epsilon[i]) / beta
  param[i] = c(alpha, beta, gamma, delta, epsilon[i], psi[i])
  eig[i] = eigen(jacobian.full(y = c(x = x[i], y = y[i]), func = eqn, parms = param[i]))$values
}

I can calculate inside the equilibrium point (vector x,y), but the function "eigen" returns "number of items to replace is not a multiple of replacement length". I guess it comes from the way I try to replace the list of parameters, I tried different ways (above is one of them) but nothing has worked. Could it be that the double function "eigen(jacobian.full(...))" does not like to be depending on indexes ?

Anybody can help ?

Solution

We can't test your code, because we don't have the jacobian.full function, but I'll guess you mean the one in the rootSolve package. When I run the code after library(rootsolve), I get these warnings:

Warning messages:
1: In param[i] <- c(alpha, beta, gamma, delta, epsilon[i], psi[i]) :
  number of items to replace is not a multiple of replacement length
2: In eig[i] <- eigen(jacobian.full(y = c(x = x[i], y = y[i]), func = eqn,  :
  number of items to replace is not a multiple of replacement length
3: In param[i] <- c(alpha, beta, gamma, delta, epsilon[i], psi[i]) :
  number of items to replace is not a multiple of replacement length
4: In eig[i] <- eigen(jacobian.full(y = c(x = x[i], y = y[i]), func = eqn,  :
  number of items to replace is not a multiple of replacement length

Those aren't coming from the eigen function, they're coming from the last two lines of your code. It's not clear what your intentions are here. param is initialized to be length 6, then you're trying to replace one element of it with a length 6 vector. Maybe the solution is simply to use

param = c(alpha, beta, gamma, delta, epsilon[i], psi[i])
eig = eigen(jacobian.full(y = c(x = x[i], y = y[i]), func = eqn, parms = param))$values

but I really can't see what you intend.