R: Решение ODE с пакетом deSolve с использованием матриц в качестве входных данных

Я пытаюсь решить обыкновенные дифференциальные уравнения в R. У меня есть матрицы в качестве начальных значений и параметров, которые я уже определил. Но когда я пытаюсь ее решить, я получаю следующую ошибку, которая не появляется, когда я ввожу отдельные значения вместо матриц.

Ошибка в (бета)% *% S: требуются числовые / комплексные матричные / векторные аргументы

Мои коды для решения оды ниже

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)

Я не понимаю, почему появляется сообщение об ошибке и нужно ли мне что-то делать по-другому, когда я использую матрицы.

Любая помощь будет принята с благодарностью. Спасибо!!!


r matrix ode
person Jessie    schedule 12.03.2017    source источник
comment
вам придется переупаковывать матрицы внутри функции градиента. Функция relist() может оказаться полезной.   -  person Ben Bolker    schedule 12.03.2017

Ответы (1)


arrow_upward
3
arrow_downward

Есть некоторые несоответствия в размерах матриц в вашем коде. Следующее должно работать (на рисунке показано, как значения параметров сходятся с итерациями):

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)

введите здесь описание изображения

person Sandipan Dey    schedule 12.03.2017
comment
@ Sandipan Dey Большое спасибо за решение. Не могли бы вы также объяснить, как матрицы относятся к функции государством? Должен ли я указывать состояние только для условий инициализации? Например, если у меня есть «S = matrix (c (S1, S2, S3), nrow = 3, ncol = 1)», то это будет «S = matrix (state [1: 3], nrow = 3, ncol = 1) в функции ?? - person Jessie; 13.03.2017
comment
Нет, потому что вы передаете параметр state, который связывается со значением трех отдельных параметров S1, S2, S3, поскольку функция не знает их по отдельности. Также «S1» означает символьную строку, а не переменную S1. - person Sandipan Dey; 13.03.2017