I would like to solve a system of differential equations with parameters varying over intervlas.
Here is my code:
# 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
The problem is that the function rigidode is working only for constant parameters. I can't figure out how to vary my parameters over interval (from 0 to 2).
thanks
@Mily comment: Yes, it is possible with t1, here the solution:
Define t1 (Intervall is not needed in my point of view).
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))
Define the ODE function:
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
Note that in the function call lsoda the argument t1 = t1 is committed to the ODE function.