R:使用矩阵作为输入，使用 deSolve 包求解 ODE

Question

我正在尝试求解 R 中的常微分方程。我有矩阵作为初始值和我已经定义的参数。但是当我尝试解决它时，出现以下错误，当我输入单个值而不是矩阵时，该错误不会出现。

Error in (beta) %*% S : requires numeric/complex matrix/vector arguments

我的解题代码如下

S = matrix(c("S1","S2"), nrow = 2, ncol=1)
I = matrix(c("I1","I2"), nrow = 2, ncol=1)
R = matrix(c("R1","R2"), nrow = 2, ncol=1)

beta=matrix(c("beta1", "beta2"), nrow = 2, ncol=1) 

MODEL <- function(time, state, parameters) {
          with(as.list(c(state, parameters)), {
          dS <- -1*(beta) %*% S %*% I
          dI <-  beta %*% S %*% I - gamma %*% I
          dR <-  gamma %*% I
          return(list(c(dS, dI, dR)))
          })
         }

init       <-c(S1=1-1e-6, S2=1-1e-6, I1=1e-6, I2=1e-6, R1=0.0, R2=0.0)
parameters <- c(beta1=1.4247, beta2=1.4247, gamma=0.14286)
times      <- seq(0, 70, by = 1)

out <- ode(y=init, times=times, func=MODEL, parms=parameters)

我不明白为什么会出现错误消息，以及我在使用矩阵时是否必须做一些不同的事情。

如有任何帮助，我们将不胜感激。谢谢!!!

Answer 1

您的代码中矩阵的维度有些不匹配。以下应该有效(图中显示了参数值如何随着迭代收敛):

MODEL <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    S = matrix(state[1:2], nrow = 2, ncol=1)
    I = matrix(state[3:4], nrow = 2, ncol=1)
    R = matrix(state[5:6], nrow = 2, ncol=1)
    beta = matrix(c(beta1, beta2), nrow = 2, ncol=1) 
    dS <- -1*(beta) %*% t(S) %*% I
    dI <-  beta %*% t(S) %*% I - gamma * I
    dR <-  gamma * I
    return(list(c(dS, dI, dR)))
  })
}

init       <- c(S1=1-1e-6, S2=1-1e-6, I1=1e-6, I2=1e-6, R1=0.0, R2=0.0)
parameters <- c(beta1=1.4247, beta2=1.4247, gamma=0.14286)
times      <- seq(0, 70, by = 1)

out <- ode(y=init, times=times, func=MODEL, parms=parameters)

library(ggplot2)
library(reshape2)
ggplot(melt(as.data.frame(as.matrix(out)), id='time'), aes(time, value, col=variable)) + 
         geom_point() + geom_line() + facet_wrap(~variable)