Решение дифференциальных уравнений с параметрами, изменяющимися с интервалами

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)) 

Определим функцию ОДУ:

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 стремится к функции ОДУ.

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 

Так что, в конце концов, я пришел к этому решению. Пожалуйста, проверьте, правильны ли все значения, которые я написал, я только что сделал ваш код!