Search code examples
rif-statementmodelingodedifferential-equations

Time varying parameter-matrix in deSolve R


I am struggling with this for so long. I have a logistic growth function where the growth parameter r is a matrix. The model is constructed in a way that I have as an output two N the N1 and N2.

I would like to be able to change the r parameter over time. When time < 50 I would like r = r1 where

r1=matrix(c(
  2,3),
  nrow=1, ncol=2

When time >= 50 I would like r=r2 where

r2=matrix(c(
  1,2),
  nrow=1, ncol=2

Here is my function. Any help is highly appreciated.

rm(list = ls())      
library(deSolve)

model <- function(time, y, params) {
  with(as.list(c(y,params)),{
    N = y[paste("N",1:2, sep = "")]
    
    dN <- r*N*(1-N/K)
    
    return(list(c(dN)))
  })
}

r=matrix(c(
  4,5),
  nrow=1, ncol=2)

K=100

params <- list(r,K)

y<- c(N1=0.1, N2=0.2)

times <- seq(0,100,1)

out <- ode(y, times, model, params)
plot(out)

I would like ideally something like this but it does not work

model <- function(time, y, params) {
  with(as.list(c(y,params)),{
    N = y[paste("N",1:2, sep = "")]
    
   r = ifelse(times < 10, matrix(c(1,3),nrow=1, ncol=2),
    ifelse(times > 10, matrix(c(1,4),nrow=1, ncol=2), matrix(c(1,2),nrow=1, ncol=2)))
     
    print(r)
    
    dN <- r*N*(1-N/K)
    
    return(list(c(dN)))
  })
}

Thank you for your time.


Solution

  • If you want to pass a matrix parameter you should pass a list of parameters and you can modify it inside the model when your time limit is exceeded (in the example below you don't even have to pass the r matrix to the model function)

    library(deSolve)
    
    model <- function(time, y, params) {
      with(as.list(c(y,params)),{
        if(time < 3) r = matrix(c(2,3), nrow = 1, ncol = 2)
        else r = matrix(c(1,3), nrow = 1, ncol = 2)
        N = y[paste("N",1:2, sep = "")]
        
        dN <- r*N*(1-N/K)
        
        return(list(c(dN)))
      })
    }
    
    y <- c(N1=0.1, N2=0.2)
    
    params <- list(r = matrix(c(0,0), nrow = 1, ncol = 2), K=100)
    times <- seq(0,10,0.1)
    
    out <- ode(y, times, model, params)
    plot(out)
    

    You can see examples of this for instance with Delay Differential Equations ?dede