1

我正在为硕士生写一个练习，我希望他们能够输入一个 ODE 系统并使用 deSolve 自动评估它。

我现在遇到的问题是根据用户输入定义要在 R 中动态评估的函数。假设我想从文本文件 ODE.txt 中读取数据，该文件包含 ODE 系统，其语法与我直接在其中写入时使用的语法相同：

  ##function (differential form)
  ODEfun <- function(time=seq(1:times_range), state=state, parms=parms) {
    
    with(as.list(c(state, parms)), {
      
      source("ODE.txt")

    })
  }
  
  ## Solve using ode (General Solver for Ordinary Differential Equations)
  out <- ode(y = init, time = seq(1:times_range), func = ODEfun, parms = parameters)

ODE.txt 文件：

      .S=-b*S*I
      .I <-  b * S * I - g * I
      .R <-                 g * I
   return(list(c(.S, .I, .R)))

（它包含我开始使用的示例系统，我已经编写了处理它的函数）

和错误：

Error in eval(ei, envir) : object 'b' not found

无论我尝试什么，错误都非常相似（我尝试了 eval(parse(, evaluate等的组合......

有谁知道如何构建这个东西？无论如何谢谢！

4

1 回答 1

3

这是一种方法。

让函数ODEfun有一个额外的参数，equations成为文本文件的内容。在函数中替换source为eval(parse()). 在ode调用之前，将文本文件读入一个变量，在下面的例子中eqs，并将该变量作为额外参数的值传递。

参数和初始值取自此RPubs 帖子。

如果答案是 parse() 你通常应该重新考虑这个问题。——Thomas Lumley R-help（2005 年 2 月））。

##function (differential form)
ODEfun <- function(time=seq(1:times_range), state=state, parms=parms, equations) {
  with(as.list(c(state, parms)), {
    
    eval(parse(text = equations))
    
  })
}

parameters <- c(
  b  = 0.004, # infectious contact rate (/person/day)
  g = 0.5    # recovery rate (/day)
)
init <- c(
  S = 999,  # number of susceptibles at time = 0
  I =   1,  # number of infectious at time = 0
  R =   0   # number of recovered (and immune) at time = 0
)
times_range <- 10

eqs <- readLines("ODE.txt")
## Solve using ode (General Solver for Ordinary Differential Equations)
out <- ode(y = init, time = 1:times_range, func = ODEfun, parms = parameters, equations = eqs)

out
#   time           S         I          R
#1     1 999.0000000   1.00000   0.000000
#2     2 963.7055522  31.79832   4.496128
#3     3 461.5686153 441.91586  96.515523
#4     4  46.1563274 569.50414 384.339532
#5     5   7.0358787 373.49828 619.465843
#6     6   2.1489403 230.12932 767.721743
#7     7   1.0390925 140.41084 858.550071
#8     8   0.6674074  85.44478 913.887809
#9     9   0.5098627  51.94497 947.545167
#10   10   0.4328912  31.56515 968.001963

编辑

这是代码 Gabor Grothendiek 的评论，稍作编辑。该函数被命名ODEfun2，其主体位于一个单独的向量中。结果identical与上述结果一致。

ODEfun2 <- function(time, state, parms) {} # empty body

fun_body <- c(
  "with(as.list(c(state, parms)), {", 
  eqs, 
  "})"
)

body(ODEfun2) <- parse(text = fun_body)
out2 <- ode(y = init, time = 1:times_range, func = ODEfun2, parms = parameters)

identical(out, out2)
#[1] TRUE

然而...

记住：

fortunes::fortune(106)

#If the answer is parse() you should usually rethink the question.
#   -- Thomas Lumley
#      R-help (February 2005)
于 2022-01-28T12:49:40.690 回答