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

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

这是（我的意思）最佳解决方案和一些解释性说明：

要使参数在函数中可用，只需使用 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

所以最后，我得出了这个解决方案。请检查我写的所有值是否正确，我刚刚让你的代码工作！

Answer 2

@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 函数。

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

solving differential equations with parameters varying over intervals

parameters

r

differential-equations