ODE, множественные корни и события, R

Question

Мне нравится решать систему связанных дифференциальных уравнений, которые включают в себя несколько порогов. Просматривая информацию R, это приводит меня к использованию ODE в сочетании с корневой функцией и функцией события.

Просматривая различные примеры, например температурную модель, стр. 14 http://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf - вставленный ниже код - я могу позволить своей модели действовать на одном пороге, то есть нахождение корня / достижение порога для одной переменной вызывает событие.

library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) c(y[1]-18, y[1]-20)

eventfunc <- function(t,y,parms) {
  y[1] <- y[1]
  y[2] <- ! y[2]  
 return(y)
}

times <- seq(from=0, to=20, by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot

В примере также показано, что разные корни могут запускать одну и ту же функцию события (y[1] - 18 и y[1]-20 оба запускают eventfunc). Мой вопрос, однако, возможно ли, чтобы разные корни запускали разные функции событий? Говорят по-разному, в зависимости от того, какой корень найден, запускается другой eventfunc? Или, в качестве альтернативы, в пределах одной и той же функции события он может выполнять различные действия в зависимости от того, какой корень найден.

Для простоты сначала я хотел посмотреть, может ли это работать, используя тот же пример, например, с помощью именования корней и работы с оператором if? На данный момент это не работает. У кого-нибудь есть опыт с этим? Если вы посмотрите на атрибуты (вне), кажется, что ODE ведет учет того, какой корень встречен, $indroot (но это после оценки.) Заранее спасибо.

# library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20
 return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  ifelse(yroot[2]==2, y[2] <- y[2], y[2] <- !y[2])
 return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot 

Answer 1

Штат системы доступен как в корневой функции, так и в функции события, поэтому его можно использовать в качестве условия для запуска события.

Для более сложных случаев, конечно, также можно отправлять события из основной функции события в другие функции для получения подробностей, то же самое можно сделать и для проверки корневого условия.

Спасибо @Bakaburg за обнаружение этого интересного вопроса без ответа.

Вот решение, которое также упрощает некоторые конструкции программирования:

      library(deSolve)
yini <- c(temp = 18, heating_on = 1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
  list(c(dy1, dy2))
}

rootfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  return(yroot)
}

eventfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
  y[2] <- if(whichroot == 2) 0 else 1
  return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
             rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)

Answer 2

Я подошел к аналогичной проблеме, используя глобальную переменную, установленную в корне или непосредственно в основной функции (полезно, если вы запускаете ее на основе преодоления порога в определенном направлении.

Затем глобальный флаг изменяет поведение функции события.

Не очень элегантно, но работает.

В вашем случае код станет таким:

      trigger <- FALSE

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20

  if (yroot[2] == 0) trigger <<- TRUE

  return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  if (trigger) y[2] <- y[2] else y[2] <- !y[2]
  return(y)
}