間隔で変化するパラメータを持つ微分方程式を解く

Question

私は、面内で変化するパラメータを持つ微分方程式の系を解きたいと思います。

# LOADING PACKAGES 
    library(deSolve) 

# DATA CREATION 
t1 <- data.frame(times=seq(from=0,to=5,by=0.1),interval=c(rep(0,10),rep(1,20),rep(2,21))) 
length(t1[which(t1$times<1),])    #10 
length(t1[which(t1$times>=1&t1$times<3),]) #20 
length(t1[which(t1$times>=3),])   #21 

t1$mueDP=c(rep(3.1,10),rep(2.6,20),rep(1.1,21)) 
t1$mueHD=c(rep(2.6,10),rep(1.7,20),rep(1.3,21)) 
t1$mueTX=c(rep(1.9,10),rep(3.3,20),rep(1.3,21)) 
t1$tau12=c(rep(5.5,10),rep(2.7,20),rep(0.7,21)) 
t1$tau13=c(rep(3.5,10),rep(1.3,20),rep(2.3,21)) 
t1$tau21=c(rep(4,10),rep(1.8,20),rep(2.8,21)) 
t1$tau23=c(rep(2.1,10),rep(2.1,20),rep(1.1,21)) 
t1$tau31=c(rep(3.9,10),rep(3.6,20),rep(1.6,21)) 
t1$tau32=c(rep(5.1,10),rep(1.4,20),rep(0.4,21)) 

t1 

# FUNCTION SOLVING THE SYSTEM 
rigidode <- function(times, y, parms) { 
with(as.list(y), { 
dert.comp_dp=-(tau12)*comp_dp+(tau21)*comp_hd-(tau13)*comp_dp+(tau31)*comp_tx-(mueDP)*comp_dp 
dert.comp_hd=-(tau21)*comp_hd+(tau12)*comp_dp-(tau23)*comp_hd+(tau32)*comp_tx-(mueHD)*comp_hd 
dert.comp_tx=-(tau31)*comp_tx+(tau13)*comp_dp-(tau32)*comp_tx+(tau23)*comp_hd-(mueTX)*comp_tx 
dert.comp_dc=(mueDP)*comp_dp+(mueHD)*comp_hd+(mueTX)*comp_tx 
list(c(dert.comp_dp, dert.comp_hd, dert.comp_tx, dert.comp_dc)) 
}) 
} 


times <- t1$times 

mueDP=t1$mueDP 
mueHD=t1$mueHD 
mueTX=t1$mueTX 
mu_attendu=t1$mu_attendu 
tau12=t1$tau12 
tau13=t1$tau13 
tau21=t1$tau21 
tau23=t1$tau23 
tau31=t1$tau31 
tau32=t1$tau32 
parms <- c("mueDP","mueHD","mueTX","mu_attendu","tau12","tau13","tau21","tau23","tau31","tau32") 
yini <- c(comp_dp = 30, comp_hd = 60,comp_tx = 10, comp_dc = 0) 

out_lsoda <- lsoda (times = times, y = yini, func = rigidode, parms = parms, rtol = 1e-9, atol = 1e-9) 
out_lsoda 

問題は、関数rigidodeは、一定のパラメータに対してのみ動作していることである。

は、ここに私のコードです。私は間隔（0から2まで）にわたって私のパラメータをどのように変えるのか分かりません。ここで

おかげ

Answer 1

@Milyコメント：

（Intervallは、私の視点では必要ありません）t1を定義します。はい、それはここでは、t1で解決可能です。関数呼び出しの引数lsodat1 = t1にODE機能にコミットしていることを

rigidode <- function(times, y, parms,t1) { 

    ## find out in which line of t1 `times` is 
    id <- min(which(times < t1$times))-1 
    parms <- t1[id,-1] 

    with(as.list(c(parms,y)), { 

    dert.comp_dp <- -(tau12)*comp_dp+(tau21)*comp_hd-(tau13)*comp_dp+(tau31)*comp_tx-(mueDP)*comp_dp 
    dert.comp_hd <- -(tau21)*comp_hd+(tau12)*comp_dp-(tau23)*comp_hd+(tau32)*comp_tx-(mueHD)*comp_hd 
    dert.comp_tx <- -(tau31)*comp_tx+(tau13)*comp_dp-(tau32)*comp_tx+(tau23)*comp_hd-(mueTX)*comp_tx 
    dert.comp_dc <- (mueDP)*comp_dp+(mueHD)*comp_hd+(mueTX)*comp_tx 

    return(list(c(dert.comp_dp, dert.comp_hd, dert.comp_tx, dert.comp_dc))) 
    }) 
} 


times <- seq(from = 0, to = 5, by = 0.1) 


yini <- c(comp_dp = 30, comp_hd = 60, comp_tx = 10, comp_dc = 0) 

parms <- t1[1,-1] 

out_lsoda <- lsoda(times = times, y = yini, func = rigidode, parms = parms, rtol = 1e-9, atol = 1e-9, t1 = t1) 
out_lsoda 

注：

t1 <- data.frame(times=seq(from=0, to=5, by=0.1)) 
t1$mueDP=c(rep(3.1,10),rep(2.6,20),rep(1.1,21)) 
t1$mueHD=c(rep(2.6,10),rep(1.7,20),rep(1.3,21)) 
t1$mueTX=c(rep(1.9,10),rep(3.3,20),rep(1.3,21)) 
t1$tau12=c(rep(5.5,10),rep(2.7,20),rep(0.7,21)) 
t1$tau13=c(rep(3.5,10),rep(1.3,20),rep(2.3,21)) 
t1$tau21=c(rep(4,10),rep(1.8,20),rep(2.8,21)) 
t1$tau23=c(rep(2.1,10),rep(2.1,20),rep(1.1,21)) 
t1$tau31=c(rep(3.9,10),rep(3.6,20),rep(1.6,21)) 
t1$tau32=c(rep(5.1,10),rep(1.4,20),rep(0.4,21)) 

はODE関数を定義します。

Answer 2

（私の意味での）最善の解決策と一部の注釈：

1. だけwith(as.list(...))機能を使用し、関数内のパラメータが利用できるようにします。

私はそれが簡単に作られたと機能に例区別を作っ：

rigidode <- function(times, y, parms) { 
    with(as.list(c(parms,y)), { 

    if(times > 1.1 & times < 3.1){  
     mueDP <- 2.6 
     mueHD <- 1.7 
     mueTX <- 3.3 
     tau12 <- 2.7 
     tau13 <- 1.3 
     tau21 <- 1.8 
     tau23 <- 2.1 
     tau31 <- 3.6 
     tau32 <- 1.4  
    } 

    if(times > 3.1){  
     mueDP <- 1.1 
     mueHD <- 1.3 
     mueTX <- 1.3 
     tau12 <- 0.7 
     tau13 <- 2.3 
     tau21 <- 2.8 
     tau23 <- 1.1 
     tau31 <- 1.6 
     tau32 <- 0.4  
    } 

    #un-comment the line below, if you want to see, if this works as expected 
    # print(c(times, mueDP, mueHD, mueTX, tau12, tau13, tau21,tau23,tau31, tau23)) 

    dert.comp_dp <- -(tau12)*comp_dp+(tau21)*comp_hd-(tau13)*comp_dp+(tau31)*comp_tx-(mueDP)*comp_dp 
    dert.comp_hd <- -(tau21)*comp_hd+(tau12)*comp_dp-(tau23)*comp_hd+(tau32)*comp_tx-(mueHD)*comp_hd 
    dert.comp_tx <- -(tau31)*comp_tx+(tau13)*comp_dp-(tau32)*comp_tx+(tau23)*comp_hd-(mueTX)*comp_tx 
    dert.comp_dc <- (mueDP)*comp_dp+(mueHD)*comp_hd+(mueTX)*comp_tx 

    return(list(c(dert.comp_dp, dert.comp_hd, dert.comp_tx, dert.comp_dc))) 
    }) 
} 

コードの残りの部分はparms回= 0の値を持っていることを、ちょうど注意し、標準です。

times <- seq(from = 0, to = 5, by = 0.1) 

yini <- c(comp_dp = 30, comp_hd = 60, comp_tx = 10, comp_dc = 0) 
parms <- c(mueDP = 3.1, mueHD = 2.6, mueTX = 1.9, tau12 = 5.5, tau13 = 3.5, 
     tau21 = 4.0, tau23 = 2.1, tau31 = 3.9, tau32 = 5.1) 

out_lsoda <- lsoda (times = times, y = yini, func = rigidode, parms = parms, rtol = 1e-9, atol = 1e-9) 
out_lsoda 

ですから、私はこの解決策になります。私が書いたすべての値が正しいかどうかを確認してください、私はちょうどあなたのコードを作った！