r - 为什么 deSolve 中的 ode 函数总是在时间 = 0 时发生事件？

Question

我正在编写一段更复杂的代码，它的行为并不像我预期的那样，经过一番修补后发现了我认为出了问题的地方:deSolve 的 ode 函数中的“事件”比我发生的次数还要多会预料到，特别是在时间 = 0 时有一个事件，即使不应该有。我在这里创建了一个最小的示例:

library(deSolve)

q <- 0

# Define the ODE system
ode_fun <- function(t, y, parms) {
  dy_dt <- -y # y' = -y
  return(list(dy_dt))
}

# Define the event function
event_fun <- function(t, y, parms) {
  print(t)
  q <<- q +1
  y <- y + 1
}

# Set up initial conditions and time points
y0 <- 0.5 # initial value of y
times <- seq(0, 5, by = 0.1) # time points to solve ODEs

# Set up the parameters for the ODE system
parms <- NULL

# Solve the ODE system with the event function
ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, times = c(2)))

# Plot the solution
plot(ode_sol, xlab = "time", ylab = "y")
print(q)

在示例中，事件应该只发生一次，所以最后 q 应该是 1，但实际上是 2! 为什么会出现这种情况？有办法解决吗？

Answer 1

这没什么好担心的。让我们调查一下会发生什么。首先，我让事件函数抛出错误以便能够看到调用堆栈:

event_fun <- function(t, y, parms) {
  print(t)
  stop("test0")
  q <<- q +1
  y + 1
}

ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
traceback()
#7: stop("test0") at #3
#6: events$func(time, state, parms, ...)
#5: Func(times[1], y)
#4: eval(Func(times[1], y), rho)
#3: checkEventFunc(Eventfunc, times, y, rho)
#2: lsoda(y, times, func, parms, ...)
#1: ode(y = y0, times = times, func = ode_fun, parms = parms, events = list(func = #event_fun, 
#       time = c(2)))

我们看到该函数是在对 checkEventFunc 的调用中被调用的。让我们看一下该函数:

deSolve:::checkEventFunc
#function (Func, times, y, rho) 
#{
#    tmp <- eval(Func(times[1], y), rho)
#    if (length(tmp) != length(y)) 
#        stop(paste("The number of values returned by events$func() (", 
#            length(tmp), ") must equal the length of the initial conditions vector #(", 
#            length(y), ")", sep = ""))
#    if (!is.vector(tmp)) 
#        stop("The event function 'events$func' must return a vector\n")
#}
#<bytecode: 0x00000173fc02a080>
#<environment: namespace:deSolve>

很明显(也从其名称来看)该函数的唯一目的是检查事件函数是否已正确定义。

现在让我们检查另一个调用:

event_fun <- function(t, y, parms) {
  print(t)
  if (t > 0) stop("test1")
  q <<- q +1
  y + 1
}

ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
traceback()
#5: stop("test1") at #3
#4: events$func(time, state, parms, ...)
#3: (function (time, state) 
#   {
#       attr(state, "names") <- Ynames
#       events$func(time, state, parms, ...)
#   })(2, 0.0676677861933153)
#2: lsoda(y, times, func, parms, ...)
#1: ode(y = y0, times = times, func = ode_fun, parms = parms, events = list(func = event_fun, 
#       time = c(2)))

我们可以看到这是对求解器实际使用的事件函数的调用。

您还可以比较有事件和没有事件的解决方案:

event_fun <- function(t, y, parms) {
  y + 1
}
ode_sol_with <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
ode_sol_without <- ode(y = y0, times = times, func = ode_fun, parms = parms)

plot(ode_sol_with, xlab = "time", ylab = "y", ylim = c(0, 1))
par(new = TRUE)
plot(ode_sol_without, xlab = "time", ylab = "y", ylim = c(0, 1), col = "red", lty = 2)

显然，在事件发生之前，两种解决方案都是相同的。