求解参数随间隔变化的微分方程

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评论：是的，这是可能的t1，这里的解决方案：

定义t1（不需要在我的观点INTERVALL）。

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功能：

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 

注意的是，在函数调用lsoda参数t1 = t1致力于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 

所以最后我来解决这个问题。请检查我写的所有值是否正确，我只是让你的代码工作！