英文:
Why does the ode function in deSolve always have an event occur at time = 0?
问题
我正在编写一段更复杂的代码，但它的行为与我的预期不符，经过一些尝试，我已经弄清楚了问题出在哪里：在deSolve的ode函数中，'events'事件的发生次数比我预期的要多一次，特别是在时间t=0时，即使不应该有事件。我在这里创建了一个最小的示例：
library(deSolve)q <- 0# 定义ODE系统ode_fun <- function(t, y, parms) {dy_dt <- -y # y' = -yreturn(list(dy_dt))}# 定义事件函数event_fun <- function(t, y, parms) {print(t)q <<- q + 1y <- y + 1}# 设置初始条件和时间点y0 <- 0.5 # y的初始值times <- seq(0, 5, by = 0.1) # 用于解ODE的时间点# 设置ODE系统的参数parms <- NULL# 使用事件函数解ODE系统ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,events = list(func = event_fun, times = c(2)))# 绘制解的图形plot(ode_sol, xlab = "时间", ylab = "y")print(q)
在这个示例中，事件应该只发生一次，所以最终q应该为1，但实际上是2！为什么会发生这种情况？有没有办法修复它？
英文:
I am writing a more complex bit of code that wasn't behaving as I expected and after some tinkering have worked out what I think is going wrong: the 'events' in the ode function of deSolve are occuring once more than I would have expected, specifically there is an event at time=0 even when there should not be. I have created a minimal example here:
library(deSolve)q <- 0# Define the ODE systemode_fun <- function(t, y, parms) {dy_dt <- -y # y' = -yreturn(list(dy_dt))}# Define the event functionevent_fun <- function(t, y, parms) {print(t)q <<- q +1y <- y + 1}# Set up initial conditions and time pointsy0 <- 0.5 # initial value of ytimes <- seq(0, 5, by = 0.1) # time points to solve ODEs# Set up the parameters for the ODE systemparms <- NULL# Solve the ODE system with the event functionode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,events = list(func = event_fun, times = c(2)))# Plot the solutionplot(ode_sol, xlab = "time", ylab = "y")print(q)
In the example, the event should only happen once, so q should be 1 at the end, but it is actually 2!
Why does this happen? Is there a way to fix it?
答案1
得分: 2
这部分代码是一个R语言的代码示例，用于模拟ODE（常微分方程）的行为，特别是与事件相关的情况。在代码中，有两个事件函数，分别命名为event_fun，它们在不同的时间点触发。
代码中的第一个event_fun函数在时间点2处抛出错误"test0"，并且会打印时间t。第二个event_fun函数在时间点2处抛出错误"test1"，也会打印时间t。这两个事件函数被用于ODE求解中，以观察它们对求解的影响。
最后，代码展示了两个ODE的解，一个包含事件（ode_sol_with），一个不包含事件（ode_sol_without），并使用plot函数可视化它们的行为。在事件发生前，两个解是相同的。
如果您需要进一步的帮助或解释，请随时提问。
英文:
This is nothing to be concerned about. Let's investigate what happens. First I let the event function throw an error to be able to see the call stack:
event_fun <- function(t, y, parms) {print(t)stop("test0")q <<- q +1y + 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)))
We see that the function was called within a call to checkEventFunc. Let's take a look at that function:
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>
It is obvious (also from its name) that the only purpose of this function is to check that the event function was defined properly.
Now let's check the other call:
event_fun <- function(t, y, parms) {print(t)if (t > 0) stop("test1")q <<- q +1y + 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)))
We can see that this is a call to the event function that is actually used by the solver.
You can also compare the solutions with and without event:
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)
Clearly, both solutions are identical before the event time.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论