【发布时间】：2017-05-01 18:18:42
【问题描述】：
我正在尝试复制现有的数学模型，以便为独立研究项目获取“控制”数据。 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)
我想我知道为什么这不起作用，但我不知道如何解决它。
有人对如何进行有任何建议吗？
【问题讨论】：
-
我认为问题在于您的模型有 5 个状态变量，因此 ODE 求解器期望 nosema 返回 5 个变量的变化率。但这不是你返回的。
pchip不会为参数生成单个值，而是为一年中的每一天生成一个值向量。通过将状态变量与这些参数相乘，可以为每个状态变量生成 365 个值。 5 * 365 = 1825。（尝试运行pchip(x, ybeta, xs)，看看你得到了什么。）我认为你的意思是pchip(x, ybeta, time)，等等。这应该会在适当的时间点产生一个参数值。 -
我认为您还需要扩展
x和相应的参数值（例如，yup），因为您只指定超过五天而不是超过一年的参数。我认为pchip无法从 5 天推断到 365 天。
标签： r parameters differential-equations