绘制微分方程组的一个解，该微分方程的系统改变R中时间的参数方面

Question

我一直在执行此代码以绘制仓本模型的解，我的目标是绘制一个改变耦合系数K的解，然后将所有图形放在一起以查看耦合的演化方式以及我在这里的尝试：

x11()
plot(times1, s1,type="l",xlim=c(0,180),col="white")
value <- seq(0, 3, by=0.1)
pe=1
for(numero in value){
  parameters <- c(N = 2, w1=0.4838171, w2=0.8842176, K=numero)
  state <- c(theta1=0, theta2=0)

  FN <- function (t, state, parameters) {
    with(as.list(c(state,parameters)),{
      dtheta1 <- 2*pi*w1 + (K/N)*(sin(theta1-theta1)+sin(theta2-theta1))
      dtheta2 <- 2*pi*w2 + (K/N)*(sin(theta1-theta2)+sin(theta2-theta2))


      list(c(dtheta1,dtheta2))
    })
  }
  times1 <- seq(1800*(value[pe]*10),1800*(value[1+pe]*10), by=1)*T
  if (pe<=29){
    pe=pe+1
  }
  library(deSolve)
  out <- ode(y = state,times = times1,func = FN,parms=parameters)
  head(out)


  s1 = (1/length(state))*(cos(out[,2])+cos(out[,3]))
  lines(times1, s1,type="l")
}

但是现在我想拥有一个随时间变化的K，它会发生变化，即像某些函数一样取K，例如K = K（t），并画出一个解决方案，使时间矢量同时改变K参数更改和随解决方案而变化的图形...但是我不知道如何解决参数随时间变化的系统，我在这里阅读了一些问题，但没有帮助，您能给我一些帮助还是如何使参数依赖于此吗？，非常感谢！

Answer 1

一般方法是定义一个用R的approxfun()插值的“强制函数”。由于您的代码中包含一些阻碍可重复性的问题，因此我将通过另一个示例进行说明。它实质上是一个逻辑增长模型，具有随时间变化的承载能力参数K：

library(deSolve)

FN <- function (t, state, parms, K_fun) {
  with(as.list(c(state, parms)), {
    K <- K_fun(t)
    dN <- r * N * (1 - N / K)
    list(c(dN))
  })
}

parms <- c(r = 1)
state <- c(N = 1)

times <- seq(0, 50)
value <- seq(10, 20, length.out = length(times))
K_fun <- approxfun(times, value, rule = 2)
out <- ode(y = state, times = times, func = FN, parms = parms, K_fun = K_fun)
plot(out)

有关此问题的更多信息，请参见软件包?events的deSolve帮助页面和以下教程：https://tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html