R 中的 deSolve ode.2D - 改变空间维度的参数

Question

我正在尝试使用 R 中 deSolve 包中的函数 ode.2D 求解二维常微分方程模型。我当前的代码是：

library(deSolve)
ModelDif <- function (time, y, pars, N, Da, dx) {
NN <- N*N
Nsp <- matrix(nrow = N, ncol = N,y[1:NN])

with (as.list(pars), {
dNsp <- Rmax * Nsp * (1- Nsp/K)
         zero <- rep(0, N)
         FluxNsp <- -Da * rbind(zero,(Nsp[2:N,] - Nsp[1:(N-1),]), zero)/dx
         dNsp    <- dNsp - (FluxNsp[2:(N+1),] - FluxNsp[1:N,])/dx
         FluxNsp <- -Da * cbind(zero,(Nsp[,2:N] - Nsp[,1:(N-1)]), zero)/dx
         dNsp    <- dNsp - (FluxNsp[,2:(N+1)] - FluxNsp[,1:N])/dx
         return(list(as.vector(dNsp)))
      })
     }
        
pars    <- c(Rmax  = 1.0,K=5)  
     R  <- 20                  
     N  <- 50                  
     dx <- R/N                   
     Da <- 0.03           
     NN <- N*N                    
     yini    <- rep(0, N*N)
     cc      <- c((NN/2):(NN/2+1)+N/2, (NN/2):(NN/2+1)-N/2)
     yini[cc] <- 1
     
     tiempo   <- seq(0, 20, by = 1)
     Final <- ode.2D(y = yini, times = tiempo, func = ModelDif, parms = pars,dimens = c(N, N), names = c("Nsp"),N = N, dx = dx, Da = Da, method = rkMethod("rk45ck"))

##Plotting
     
P1 <- subset(Final, select = "Nsp", arr = TRUE)
for(i in 1:length(tiempo)){
TT<-tiempo[i]
image(as.matrix(P1[,,i]), xlab = "Lat", ylab = "Lon")
mtext(paste("tiempo = ", TT),side=3)
Sys.sleep(0.5)}

但是，现在我需要将主方程中的固定参数 K 替换为从与用于求解模型的网格大小相同的矩阵中的各个单元格获得的值。

例如，将 K 的固定值替换为 K<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N) 并在模型的每次评估中将 K 替换为相应单元格中的值的等式。

我已经在站点中搜索解决方案，但 none 有效。我会很感激任何建议。谢谢

Answer 1

如果您使用 K<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N)，则需要将所有零替换为正数，因为您将其用作分母。

在下面的模型中，我将 0 替换为 1。您可以将其替换为您选择的其他数字。

ModelDif <- function (time, y, pars, N, Da, dx) {
  NN <- N*N
  Nsp <- matrix(nrow = N, ncol = N,y[1:NN])
  KK<-matrix(rbinom(N*N,10,0.3),ncol=N,nrow=N)
  for (i in 1:NN) {
    if (round(KK[i],1) == 0) {KK[i] <- 1}  ###  replace all zeros with 1
  }
  NspK <- Nsp/KK
  
  with (as.list(pars), {
    dNsp <- Rmax * Nsp * (1- NspK)
    zero <- rep(0, N)
    FluxNsp <- -Da * rbind(zero,(Nsp[2:N,] - Nsp[1:(N-1),]), zero)/dx
    dNsp    <- dNsp - (FluxNsp[2:(N+1),] - FluxNsp[1:N,])/dx
    FluxNsp <- -Da * cbind(zero,(Nsp[,2:N] - Nsp[,1:(N-1)]), zero)/dx
    dNsp    <- dNsp - (FluxNsp[,2:(N+1)] - FluxNsp[,1:N])/dx
    return(list(as.vector(dNsp)))
  })
}

此外，运行需要一段时间，并给出以下警告： rk(y, times, func, parms, method = method, ...) 中的警告：时间步数 104651 在 t = 16.6229

时超过 maxsteps

您可能需要进一步调整一些参数。

R 中的 deSolve ode.2D - 改变空间维度的参数

deSolve ode.2D in R - changing parameter in spatial dimensions

r

matrix

ode

cran