1

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

我想将三个微分方程与 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 链接

4

1 回答 1

0

代码像这样完美运行

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)
于 2021-07-10T10:58:57.757 回答