4

私は、複数のしきい値を含む連成微分方程式系を解くのが好きです。R の情報を調べると、ODE をルート関数およびイベント関数と組み合わせて使用​​することがわかります。

さまざまな例、つまり温度モデル、14 ページhttp://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf -- 以下にコードを貼り付けます --、モデルを機能させることができますつまり、ルートを見つけたり、1 つの変数のしきい値に達したりすると、イベントがトリガーされます。

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 をトリガーします)。私の質問は、しかし、異なるルートが異なるイベント関数をトリガーすることも可能ですか? 別の言い方をすると、見つかったルートに応じて、別のイベント関数がトリガーされますか? または、同じ eventfunct 内で、見つかったルートに応じて異なるアクションを実行できます。

簡単にするために、たとえばルートに名前を付けたり、if ステートメントを操作したりするなど、同じ例を使用してこれが機能するかどうかを最初に確認したかったのですか? 現時点では機能しません。誰もこれを経験していますか?attributes(out) を見ると、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 
4

2 に答える 2