deSolve R中的时变参数矩阵

Question

我为这件事挣扎了这么久。我有一个logistic增长函数，其中增长参数r是一个矩阵。该模型的构造方式，我有一个输出，两个N，N1和N2。

我希望能随着时间的推移改变r参数。当时间< 50时，我想要r= r1

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

当时间>= 50，我想r=r2在哪里

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

这是我的功能。任何帮助都是非常感谢的。

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)

我想要这种理想的东西，但它不起作用

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

谢谢您抽时间见我。

Answer 1

如果您想传递一个矩阵参数，您应该传递一个参数列表，并且在超过时间限制时，您可以在模型内修改它(在下面的示例中，您甚至不必将r矩阵传递给模型函数)

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)

您可以看到这方面的例子，例如延迟微分方程?dede。

Answer 2

这里是一种通用方法，它使用approx函数的扩展版本。还请注意模型函数和参数值的附加图的进一步简化。

编辑根据刘易斯·卡特的建议在t=3上修改参数，这样就可以看到效果。

library(simecol) # contains approxTime, a vector version of approx

model <- function(time, N, params) {
    r <- approxTime(params$signal, time, rule = 2, f=0, method="constant")[-1]
    K <- params$K

    dN <- r*N*(1-N/K)
    
    return(list(c(dN), r))
}

signal <- matrix(
  # time, r[1, 2], 
  c(  0, 2, 3,
      3, 1, 2,
    100, 1, 2), ncol=3, byrow=TRUE
) 

## test of the interpolation
approxTime(signal, c(1, 2.9, 3, 100), rule = 2, f=0, method="constant")

params <- list(signal = signal, K = 100)

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

times <- seq(0, 10, 0.1)

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

​

​

对于示例中的少量状态变量，使用approxfun与包stats分离的信号看起来不那么通用，但速度可能更快。

作为进一步的改进，人们可以考虑用一个更平稳的过渡来代替“艰难”的过渡。然后，可以直接将其表述为不需要approx、approxfun或approxTime的函数。

编辑2:

Package simecol导入deSolve，我们只需要从它获得一个小的函数。因此，不必加载simecol，也可以在代码中显式地包含approxTime函数。从数据帧到矩阵的转换提高了性能，但是在这种情况下，矩阵仍然是首选的。

approxTime <- function(x, xout, ...) {
  if (is.data.frame(x)) {x <- as.matrix(x); wasdf <- TRUE} else wasdf <- FALSE
  if (!is.matrix(x)) stop("x must be a matrix or data frame")
  m <- ncol(x)
  y <- matrix(0, nrow=length(xout), ncol=m)
  y[,1] <- xout
  for (i in 2:m) {
    y[,i] <- as.vector(approx(x[,1], x[,i], xout, ...)$y)
  }
  if (wasdf) y <- as.data.frame(y)
  names(y) <- dimnames(x)[[2]]
  y
}