Почему собственная функция возвращает «количество заменяемых элементов не кратно длине замены»?

Question

Я использую для вопроса здесь классический пример модели Лотки-Вольтерры с хищничеством (2 ОДУ, 6 параметров). Мне нужно вычислить точку равновесия, в которой я знаю аналитическое выражение, и собственные значения матрицы Якоби в этой точке. Мне нужно сделать это для большого количества местоположений (здесь я сделаю только 2), которые отличаются в этом примере значением 2 параметров, epsilon и psi (2 значения параметров для каждого для 2 местоположений).

Я создал цикл (потому что опять размер эпсилон и пси будет намного больше). Вот мой код:

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
}

Я могу вычислить внутри точки равновесия (вектор x, y), но функция «eigen» возвращает «количество заменяемых элементов не кратно длине замены». Я думаю, это происходит из-за того, что я пытаюсь заменить список параметров, я пробовал разные способы (выше один из них), но ничего не помогло. Может ли быть так, что двойная функция "eigen (jacobian.full (...))" не любит зависеть от индексов?

Кто-нибудь может помочь?

Answer 1

Мы не можем протестировать ваш код, потому что у нас нет функции jacobian.full, но я думаю, вы имеете в виду ту, что в пакете rootSolve.Когда я запускаю код после library(rootsolve), я получаю следующие предупреждения:

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

Они приходят не из функции eigen, а из двух последних строк вашего кода.Не ясно, каковы ваши намерения здесь.param инициализируется длиной 6, затем вы пытаетесь заменить один его элемент вектором длины 6.Может быть, решение состоит в том, чтобы просто использовать

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

, но я действительно не понимаю, что вы собираетесь.