Блог

Я пытаюсь решить модель обыкновенных дифференциальных уравнений в 2D, используя функцию ode.2D в пакете отменяется в R. Мой текущий код:

 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<-матрицу (rbinom(N * N,10,0.3), ncol=N,nrow =N) и при каждой оценке модели заменяйте K в уравнении значением в соответствующей ячейке.

У меня есть поиск решения на сайте, но ни одно из них не сработало. Я буду признателен за любой совет. Спасибо

Если вы используете 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, метод = метод, …) :
Количество временных шагов 104651 превысило maxsteps при t = 16.6229

Возможно, вам потребуется дополнительно настроить некоторые параметры.