我有以下3个主要代码片段，首先绘制降雨量，然后绘制降雨量对猎物(资源)增长率的影响，然后绘制使用恒定增长率的消费者-资源(食草动物-植物)动态。我的目标是在消费者资源模型中实现动态增长率(通过方程dg)。我的最终目标是一个消费者-资源人口随时间动态增长率的动态图表。
第一段代码(绘制降雨量):
### Rainfall plot ###
t = seq(0, 50, 1) # time period
avgrain = 117.4 # average rainfall
A = 100
w = 0.6
phi = 0.1
rain = avgrain + (A*sin((w*t)+phi)) # rainfall function
plot(t, rain, type="l", xlab="time", ylab="Rainfall") # rainfall plot
第二段代码(绘制降雨量对资源增长率的影响):
### Rainfall's effect on growth rate, g ###
ropt1 = 117.4 # optimal rainfall for resource growth
s1 = 120 # standard deviation for resource growth rate as a function of rainfall
dg = exp(-(rain - ropt1)^2/(s1^2)) # rain's effect on plant growth rate
plot(t, dg, type="l", xlab="time", ylab="Plant growth rate as a function of rainfall")
第三段代码(绘制消费者资源动态图——这是我试图实现上面创建的动态增长率dg，而不是静态增长率g的地方):
### Consumer-resource model as a function of time ###
library(deSolve)
states <- c(r=1, # resource (plant population) state variable
c=1) # consumer (herbivore population) state variable
parameters <- c(g=1, # resource growth rate )
K=25, # resource carrying capacity
a=0.5, # consumer attack rate (between 0-1)
h=1, # consumer handling time
e=0.9, # consumer conversion efficiency
m=0.5) # consumer mortality rate
function1 <- function(times1, states, parameters) {
with(as.list(c(states, parameters)),{
# rate of change of state variables
dr = g*r*(1-(r/K))-((c*a*r)/(1+(a*h*r)))
dc = ((c*e*a*r)/(1+(a*h*r)))- c*m
# return rate of change
list(c(dr, dc))
})
}
times1 <- seq(0, 100, by = 1)
out1 <- ode(y = states, times = times1, func = function1, parms = parameters, method="ode45")
plot(out1) # plot state variable change across time
因此，本质上，在每个时间步长，我希望消费者资源动态根据该时间步长的增长率进行更新。这可能吗?如果有，怎么做?提前感谢您的友好回复。
可以直接将降雨方程包含在model函数中。这里t总是实际的时间步长。它们在图中也作为"辅助变量"出现。如果我们将它们作为返回值的第三和第四个元素加在导数向量之后。
若dg为无量纲值，可与增长率相乘，若为增长率，则将g替换为dg。请检查这个!
来自降雨子模型的附加参数可以保持全局，也可以在参数向量中移动。我还重命名了一些标识符(以1结尾的标识符)，并将求解器更改为lsoda而不是ode45。
library(deSolve)
states <- c(r=1, # resource (plant population) state variable
c=1) # consumer (herbivore population) state variable
parameters <- c(g=1, # resource growth rate )
K=25, # resource carrying capacity
a=0.5, # consumer attack rate (between 0-1)
h=1, # consumer handling time
e=0.9, # consumer conversion efficiency
m=0.5, # consumer mortality rate
s1=120,
avgrain = 117.4, # average rainfall
A = 100,
w = 0.6,
phi = 0.1,
ropt1 = 117.4, # optimal rainfall for resource growth
s1 = 120 # sd for resource growth rate ...
)
# note that "t" is only one time point from times
model <- function(t, states, parameters) {
with(as.list(c(states, parameters)), {
# rainfall time series function
rain <- avgrain + (A*sin((w*t)+phi)) # rainfall function
# functional response equation
dg <- exp(-(rain - ropt1)^2/(s1^2)) # rain's effect on plant growth rate
# rate of change of state variables
dr <- dg * g*r*(1-(r/K)) - ((c*a*r)/(1+(a*h*r)))
dc <- ((c*e*a*r)/(1+(a*h*r)))- c*m
# return rate of change
list(c(dr, dc), rain=rain, dg=dg)
})
}
times <- seq(0, 100, by = 1)
## lsoda is faster and more robust that ode45
out <- ode(y = states, times = times, func = model, parms = parameters, method="lsoda")
plot(out) # plot states and auxiliary variables across time