终止涉及根函数和事件函数的 ODE 求解器

Question

出于演示目的，我从文档中获取了这个示例。在以下示例中，当 y 的值达到 0.1 时，将添加一个随机值。如果 y 值大于 0.8，我想终止求解器。

一种可能的解决方案是在 eventfun 中生成一个随机值，使 y 始终小于 0.8。还有其他可能的解决方案来终止求解器吗？这对我的复杂模型很有帮助。

## =======================================================================
## Example 3:
##   using lsodar to trigger an event
## =======================================================================

## a state variable is decaying at a first-order rate. 
## when it reaches the value 0.1, a random amount is added.
library("deSolve")

derivfun <- function (t,y,parms)
  list (-0.05 * y)

rootfun <- function (t,y,parms)
  return(y - 0.1) 

eventfun <- function(t,y,parms)
  return(y + runif(1))  


yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, events = list(func=eventfun, root = TRUE))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

# }

Answer 1

以下是您想要的吗？感谢您提供清晰的示例。

library("deSolve")

derivfun <- function (t,y,parms)
  list (-0.05 * y)

rootfun <- function (t,y,parms)
  return(c(y - 0.1, y - 0.8))

eventfun <- function(t,y,parms)
  return(y + runif(1))  


yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, 
              events = list(func=eventfun, root = TRUE, terminalroot = 2))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

这里是 rootfun 和 eventfun 的另一个改进：

library("deSolve")

terminate <- 0.8
eps       <- 1e-6

derivfun <- function (t, y, parms)
  list (-0.05 * y)

rootfun <- function (t, y, parms)
  return(c(y - 0.1, y - terminate))

eventfun <- function(t, y, parms)
  return(min(y + runif(1), terminate + eps))  

yini <- 0.5
times <- 0:400

out <- lsodar(func=derivfun, y = yini, times=times, 
              rootfunc = rootfun, 
              events = list(func = eventfun, root = TRUE, terminalroot = 2))

plot(out, type = "l", lwd = 2, main = "lsodar with event")

终止涉及根函数和事件函数的 ODE 求解器

Terminate ODE solver involving root function and event function

r

ode