如何解决R中具有时间相关参数的常微分方程？

Question

我正在尝试复制现有的数学模型，以便为独立研究项目获得“控制”数据。 Link到文章。

我想使用deSolve软件包求解R中包含五个常微分方程的系统。但是，大多数参数是季节性的。先前的研究人员使用了pracma软件包中的“ pchip”函数来创建参数函数：

pchip是一种“形状保持”的分段三次Hermite多项式方法，它试图确定斜率，以使函数值不会超过数据值。 pchipfun是pchip的包装器，并返回一个函数。 pchip和pchipfun返回的函数都被向量化。

xi和yi必须是长度相等或大于3的向量（可能进行三次插值），并且xi必须进行排序。 pchip可以应用于[min（xi），max（xi）]以外的点，但是在此间隔之外的结果没有太大意义。

graphed, 'beta' parameter over five years

这是我的代码（只是试图获得一年的输出）：

x <- c(0, 1, 2, 3, 4)
ybeta <- c(500, 1500, 500, 0, 500)
yk <- c(8000, 12, 0, 8000, 8000)
yrec <- c(0.25, 0.25, 0.25, 0, 0.25)
yfb <- c(1.5, 1.5, 1.5, 1.5, 1.5)
yno <- c(0, 0, 0, 0.00649, 0)
yni <- c(0, 0, 0, 0.00649, 0)
ypo <- c(0.08511, 0.08511, 0.08511, 0, 0.08511)
ypi <- c(0.16936, 0.16936, 0.16936, 0, 0.16936)
yalpha <- c(0.55, 0.12, 0.24, 0, 0.55)
yW <- 0.1
ydep <- c(0.2061 * yW, 0.2835 * yW, 0.2527 * yW, yW, 0.2061 * yW)
ydec <- c(0.006470, 0.023300, 0.015683, 0, 0.006470)
yup <- c(0.15, 0.15, 0.15, 0.15, 0.15)

xs <- seq(0, 1, len = 365)


nosema <- function(time, state, parameters) {

    with(as.list(c(state, parameters)), {

        H <- Ho + Hi
        Fr <- Fo + Fi
        Z <- H + Fr
        dHo <- pchip(x, ybeta, xs)* (Z^n) / (pchip(x, yk, xs)^n + Z^n) - pchip(x, yrec, xs) * Ho + pchip(x, yfb, xs) * Fr/Z * Fo - pchip(x, yno, xs) * Ho - pchip(x, yalpha, xs) * Ho * E/(sr + E)
        dHi <- -pchip(x, yrec, xs) * Hi + pchip(x, yfb, xs) * Fr/Z * Fi - pchip(x, yni, xs) * Hi + pchip(x, yalpha, xs) * Ho * E/(sr + E)
        dFo <- pchip(x, yrec, xs) * Ho - pchip(x, yfb, xs) * Fr/Z * Fo - pchip(x, ypo, xs) * Fo
        dFi <- pchip(x, yrec, xs) * Hi - pchip(x, yfb, xs) * Fr/Z * Fi - pchip(x, ypi, xs) * Fi
        dE <- pchip(x, ydep, xs) * Hi - pchip(x, ydec, xs) * E - pchip(x, yup, xs) * Ho * E / (sr + E)

        return(list(c(dHo, dHi, dFo, dFi, dE)))
    })
}


init <- c(Ho = 10^4, Hi = 0, Fo = 10000, Fi = 0, E = 0)

parameters <- c(n = 2, 
        sr = 10000)

out <- ode(y = init, times = seq(0, 365, len = 365), func = nosema, parms = parameters)

out <- as.data.frame(out)

out$time <- NULL

head(out)

但是，当我运行上面的代码时...

由func（）返回的导数数（1825）必须等于初始条件向量（5）的长度

[我想我对这为什么不起作用有一个暗示，但是我不知道如何解决它。

有人对如何进行有任何建议吗？

Answer 1

library(deSolve)
?forcings

请参见第72ff页的http://desolve.r-forge.r-project.org/user2014/tutorial.pdf