R:Solving ODE with deSolve package using matrices as input

Question

I'm trying to solve ordinary differential equations in R. I have matrices as initial values and parameters which I have already defined. But when I try to solve it I'm getting the below error which is not appearing when I input single values instead of matrices.

Error in (beta) %*% S : requires numeric/complex matrix/vector arguments

My codes for solving ode is below

S = matrix(c("S1","S2"), nrow = 2, ncol=1)
I = matrix(c("I1","I2"), nrow = 2, ncol=1)
R = matrix(c("R1","R2"), nrow = 2, ncol=1)

beta=matrix(c("beta1", "beta2"), nrow = 2, ncol=1) 

MODEL <- function(time, state, parameters) {
          with(as.list(c(state, parameters)), {
          dS <- -1*(beta) %*% S %*% I
          dI <-  beta %*% S %*% I - gamma %*% I
          dR <-  gamma %*% I
          return(list(c(dS, dI, dR)))
          })
         }

init       <-c(S1=1-1e-6, S2=1-1e-6, I1=1e-6, I2=1e-6, R1=0.0, R2=0.0)
parameters <- c(beta1=1.4247, beta2=1.4247, gamma=0.14286)
times      <- seq(0, 70, by = 1)

out <- ode(y=init, times=times, func=MODEL, parms=parameters)

I do not understand why the error message appears and whether I have to do something differently when I use matrices.

Any help would be greatly appreciated. Thanks!!!

you're going to have to "repack" your matrices inside your gradient function. The relist() function might be useful. — Ben Bolker

Sandipan Dey Sandipan Dey · Accepted Answer · 2017-03-12T08:51:30

There are some mismatched in the dimensions of the matrices in your code. The following should work (the figure shows how the parameter values converge with iterations):

MODEL <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    S = matrix(state[1:2], nrow = 2, ncol=1)
    I = matrix(state[3:4], nrow = 2, ncol=1)
    R = matrix(state[5:6], nrow = 2, ncol=1)
    beta = matrix(c(beta1, beta2), nrow = 2, ncol=1) 
    dS <- -1*(beta) %*% t(S) %*% I
    dI <-  beta %*% t(S) %*% I - gamma * I
    dR <-  gamma * I
    return(list(c(dS, dI, dR)))
  })
}

init       <- c(S1=1-1e-6, S2=1-1e-6, I1=1e-6, I2=1e-6, R1=0.0, R2=0.0)
parameters <- c(beta1=1.4247, beta2=1.4247, gamma=0.14286)
times      <- seq(0, 70, by = 1)

out <- ode(y=init, times=times, func=MODEL, parms=parameters)

library(ggplot2)
library(reshape2)
ggplot(melt(as.data.frame(as.matrix(out)), id='time'), aes(time, value, col=variable)) + 
         geom_point() + geom_line() + facet_wrap(~variable)

R:Solving ODE with deSolve package using matrices as input

1 Answers