eval 中的错误（expr、envir、enclos）：找不到对象“B”答案

【问题标题】：Error in eval(expr, envir, enclos) : object 'B' not foundeval 中的错误（expr、envir、enclos）：找不到对象“B”
【发布时间】：2017-04-02 15:20:07
【问题描述】：

我被这段代码困住了将近两个月，任何帮助都将不胜感激。

我想将三个微分方程与 R 中的 deSolve 包集成。这是我的代码

library(deSolve)
library(ggplot2)


### Parameters 

D = 0.1  
S0= 6 
c = 2.3 * 10 ^-5 
a = c (0.25, 0.225, 0.2, 0.175, 0.15) # algae maximum growth rate 
H = 1 # algae conversion efficiency 
phi = 7.5 * 10^-8  
beta = 100  
epsilon = 10^-3

M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon),    (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
          nrow=5,
          ncol=5,
          byrow=TRUE)

 M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
          nrow=4,
          ncol=4,
          byrow=TRUE)



 A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
      nrow=5, 
      ncol=4,
      byrow=TRUE)

  ## time sequence
 time <- seq(0,1000, by = 1)

# parameters: a named vector
parameters <- c(D = 0.1,
            c = 2.3, 
            H = 1,
            a = c (0.25, 0.225, 0.2, 0.175, 0.15),
            S0= 30, 
            c = 2.3 * 10 ^-5,  
            H = 1, 
            phi = 7.5 * 10^-8,  
            beta = 100,
            epsilon = 10^-3,
            M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon), (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
                          nrow=5,
                          ncol=5,
                          byrow=TRUE),
            M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
                          nrow=4,
                          ncol=4,
                          byrow=TRUE),
            A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
                        nrow=5, 
                        ncol=4,
                        byrow=TRUE)) 


    nutrients <- function(t, state, parameters){
 with(as.list(c(state, parameters)),{
 g= a*S / (H + S)
 dS= D*(S0 - S) - c*sum(g,B)
 dB = M_B %*% (g * B) - (phi * (A %*% P)) * B - D*B
 dP= (M_P * beta) %*% (phi*(t(A)%*%B)*P) - (phi*(t(A)%*%B)*P) - D*P   
 return(list(c(dS,dB,dP)))
 })
 }

  out <- ode(y = c(S=30, B=c(10000,0,0,0,0), P=c(100,0,0,0)), times = time,    func = nutrients, parms = parameters)

但是，自从我收到此错误以来，我还没有成功：

eval(expr, envir, enclos) 中的错误：找不到对象“B”

你知道我做错了什么吗？

更新

经过一段时间的尝试，我找到了问题的答案。稍后我将发布一个 github 链接，其中包含解决方案和图表

【问题讨论】：

B 和 P 应该是数字还是矩阵？还是它们是 5 维和 4 维的向量？但是像(phi*t(A)*B)*P 这样的构造就没有意义了，...
@LutzL B 和 P 应该分别是 5,1 和 4,1 维度的矩阵。你很可能有一个观点，这个“（phi * t（A）* B）* P”没有意义，但我不明白。我唯一要做的就是用 R 写出一篇已经发表的论文中的微分方程......为什么没有意义:)？
我得到的错误是“错误：找不到对象'a'”，我认为你想让参数成为一个列表而不是一个原子向量。看起来S 没有在 g 表达式可以找到它的任何地方定义。如果您想获得信息性错误消息，需要从“空白石板”开始清除您的工作区。
那为什么你的初始值只是数字？ --我不知道R，所以请帮我解释一下%*%应该做什么操作。 -- 你能以伪数学的方式写下原始方程或以其他方式记录它们吗？
请给我们引用您想要关注的论文。

标签： r ode differential-equations desolve

【解决方案1】：

代码像这样完美运行

library(deSolve);library(ggplot2);library(reshape2) # load packages deSolve is for solving ODE

##############################################################################################################
#
# model parameters
#
##############################################################################################################

a = c(0.25, 0.225, 0.2, 0.175, 0.15) #alga growth rate for each one of the five types
epsilon = 10^-3                      #mutation rate

M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon), (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
              nrow=5,
              ncol=5,
              byrow=TRUE)           #mutation-transition matrix for the host the sum of each column should be one


M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
              nrow=4,
              ncol=4,
              byrow=TRUE)            #mutation-transition matrix for the virus the sum of each column should be one


A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
          nrow=5, 
          ncol=4,                   #infection matrix where 1 means infection and 0 resistance
          byrow=TRUE)


# here I put all the parameters into a vector 
parameters <- c(D = 0.1,
                epsilon = 10^-3,
                a = c(0.25, 0.225, 0.2, 0.175, 0.15),
                S0 = 10, # inflow rate 
                c = 2.3 * 10 ^-5, # algae conversion coefficiency 
                H = 1, # algae conversion efficiency 
                phi = 7.5 * 10^-8,  #  virus adsorption rate 
                beta = 200, # virus burst size 
                M_B =  matrix(c(1-epsilon, epsilon/2, 0,0,0,epsilon, (1-epsilon), (epsilon/2), 0, 0, 0, epsilon/2, (1-epsilon), epsilon/2, 0, 0, 0 , epsilon/2, (1-epsilon), epsilon,0,0,0, epsilon/2, 1-epsilon),
                              nrow=5,
                              ncol=5,
                              byrow=TRUE),
                M_P =  matrix(c(1-epsilon, epsilon/2,0,0,epsilon, (1-epsilon),(epsilon/2), 0, 0, epsilon/2, (1-epsilon), epsilon, 0,0,  epsilon/2, (1-epsilon)),
                              nrow=4,
                              ncol=4,
                              byrow=TRUE),
                A= matrix(c(1,1,1,1,0,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0), 
                          nrow=5, 
                          ncol=4,
                          byrow=TRUE)) 


##############################################################################################################
#
# time
#
##############################################################################################################

time <- seq(0,300, by = 0.01)

##############################################################################################################
#
# initial values
#
##############################################################################################################

yini= c(S=30, B=c(100,0,0,0,0), P=c(1,0,0,0))

##############################################################################################################
#
# function
#
##############################################################################################################


nutrients <- function(t, yini, parameters) {
  with(as.list(c(yini, parameters)), {
    g = ((a*S) / (H + S))
    B = yini[paste("B",1:5, sep = "")]   #### like this you convert B from a vector to an array and this is what the program needs 
    P = yini[paste("P",1:4, sep = "")]
    
    dS = (D*(S0 - S) - c*sum(g,B))
    dB = M_B %*% (g * B) - (phi * (A %*% P)) * B - D*B 
    dP = (M_P * beta) %*% (phi*(t(A)%*%B)*P) - (phi*(t(A)%*%B)*P) - D*P   
    
    
    return(list(c(dS,dB,dP)))  
  })
}

##############################################################################################################
#
# solve the model
#
##############################################################################################################

out <- ode(y = yini, times = time, func = nutrients, parms = parameters, method="lsoda",atol=10^-15,rtol=10^-15)

【讨论】：