solving differential equations with parameters varying over intervals

Question

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

You should solve it piecewise, jumps in the ODE function play havoc with adaptive step size algorithms like lsoda. Thus solve from start to first jump point, change the constants, take the last state as the initial state and solve from first to second jump point etc. — Lutz Lehmann

J_F J_F · Accepted Answer · 2016-09-28T15:50:11

Here the (in my meaning) best solution and some explanatory notes:

To make parameters available in the function, just use the with(as.list(...)) function.

I made it easy and made a distinction of cases in the function:

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

The rest of the code is standard, just note, that the parms have the values of the times = 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

So in the end, I come to this solution. Please check if all the values I wrote are right, I just made your code work!

solving differential equations with parameters varying over intervals

2 Answers