Написание дискретной системы диффузионно-реакционной модели в R

Question

У меня есть диффузионно-реакционная (самоорганизация растительности) дискретная система уравнений с заданными параметрами. Я установил начальные и граничные условия. Я считаю, что циклы for написаны правильно, но у меня ошибка. Пожалуйста, мне нужна помощь, что может быть не так.

cb <- 10 # water to biomass density conversion effeciency in [g/m2mm]
db <- 0.25 #biomass mortality rate in [1/d]
kb <- 5 #half sat. constant of specific growth and uptake 
gb <- 0.05 #max. specific water uptake
Db <- 0.1 #diffusion rate of vegetation in [m2/d]
gamma <- 0.2 #max. infiltration coeff. in [1/d]
rw <- 0.2 # soil water losses(evaporation/percolation) in [1/d]
Dw <- 0.1 #diffusion coeff. of soil water in [m2/d]
Dh <- 10 #advection coeff. for downslope water flow in [m/d]
R <- 1 # rainfall intensity in [mm/d]
k <- 5 #saturation constant for water infiltration 
W0 <- 0.2 # water infiltration rate without plants
L <- 10 # length in [m]
Tf <- 10 #duration in days
dt <- 1 #time step in days
dx <-1  #space step in [m]
Nt <- Tf/dt
Nx <- L/dx
t <- numeric(Nt)
B <- matrix(nrow = Nt, ncol = Nx) #plant biomass density
W <- matrix(nrow = Nt, ncol = Nx) #soil water
H <- matrix(nrow = Nt, ncol = Nx) #surface water
U <- matrix(nrow = Nt, ncol = Nx) #Water uptake
I <- matrix(nrow = Nt, ncol = Nx) #Soil water infiltration
B[1,]<-0
B[1,1]<-1
W[1,]<-0
W[1,3]<-1
H[1,]<-0
H[1,]<-1
U[1,]<-0
U[1,2]<-1
I[1,]<-0
I[,2]<-1

for(i in 1:Nt){
  t[i] <- i*dt
  for(j in 2:Nx-1){  

     Bleft <- B[i, j-1]
     Bright <- B[i, j+1]  
     if (j==1){
     Bleft<- B[i,Nx]
     }
     if (j==Nx){
     Bright  <- B[i,1]
     }    
     Wleft <- W[i, j-1]
     Wright <- W[i, j+1]    
     if (j==1){
     Wleft <- W[i,Nx]
     }
     if (j==Nx){
     Wright <- W[i,1]
     }  
     Hleft <-  H[i, j-1]
     Hright <- H[i, j+1] 
     if (j==1){
     Hleft <- H[i,Nx]
     }
     if (j==Nx){
     Hright <- H[i,1]
     }   
     U[i,j] <- gb*B[i,j]*W[i,j]/(W[i,j] + kb)

     I[i,j] <- gamma*H[i,j]*(B[i,j] + k*W0)/(B[i,j] + k)

     B[i+1,j] <- floor(B[i,j] + dt*((cb*U[i,j]) - (db*B[i,j]) + 
     (Db/(dx^2))*(Bleft - 2*B[i,j] + Bright)))

     W[i+1,j] <- floor(W[i,j] + dt*((I[i,j]  - U[i,j] - rw * 
     W[i,j]) + (Dw/(dx^2))*(Wleft - 2*W[i,j] + Wright)))

     H[i+1,j] <- floor(H[i,j] + dt*((R - I[i,j]) + (Dh/(dx^2))* 
     (Hleft - 2*H[i,j] + Hright)))

   }
  }

### Ошибка в `[<-` (` <em>tmp `, i + 1, j, значение = NA_real_):
нижний индекс вне границ

Answer 1

Чтобы дать вам отправную точку?

cb <- 10 # water to biomass density conversion effeciency in [g/m2mm]
db <- 0.25 #biomass mortality rate in [1/d]
kb <- 5 #half sat. constant of specific growth and uptake 
gb <- 0.05 #max. specific water uptake
Db <- 0.1 #diffusion rate of vegetation in [m2/d]
gamma <- 0.2 #max. infiltration coeff. in [1/d]
rw <- 0.2 # soil water losses(evaporation/percolation) in [1/d]
Dw <- 0.1 #diffusion coeff. of soil water in [m2/d]
Dh <- 10 #advection coeff. for downslope water flow in [m/d]
R <- 1 # rainfall intensity in [mm/d]
k <- 5 #saturation constant for water infiltration 
W0 <- 0.2 # water infiltration rate without plants
L <- 10 # length in [m]
Tf <- 10 #duration in days
dt <- 1 #time step in days
dx <-1  #space step in [m]
Nt <- Tf/dt
Nx <- L/dx
t <- numeric(Nt)

УСТАНАВЛИВАЕТСЯ НА НОЛЬ ИНИЦИАЛИЗИРОВАННЫЙ ВМЕСТО НА ПЕРХАПЕ НА

B <- matrix(0, nrow = Nt, ncol = Nx) #plant biomass density
W <- matrix(0, nrow = Nt, ncol = Nx) #soil water
H <- matrix(0, nrow = Nt, ncol = Nx) #surface water
U <- matrix(0, nrow = Nt, ncol = Nx) #Water uptake
I <- matrix(0, nrow = Nt, ncol = Nx) #Soil water infiltration

B[1,]<-0
B[1,1]<-1
W[1,]<-0
W[1,3]<-1
H[1,]<-0
H[1,]<-1
U[1,]<-0
U[1,2]<-1
I[1,]<-0
I[,2]<-1

А потом я почистил модель и добавил круглые скобки - вы можете выяснить, что здесь происходит не так, я уверен?

for(i in 1:Nt){
  t[i] <- i*dt
  for(j in 2:(Nx-1)){  # NOTICE THE PARENTHESES

    Bleft <- B[i, j-1]; Bright <- B[i, j+1]  
    Wleft <- W[i, j-1]; Wright <- W[i, j+1]      
    Hleft <- H[i, j-1]; Hright <- H[i, j+1]     

    if (j==1){
      Bleft <- B[i,Nx]
      Wleft <- W[i,Nx]
      Hleft <- H[i,Nx]
    } 
    if (j==Nx){
      Bright <- B[i,1]
      Wright <- W[i,1]
      Hright <- H[i,1]
    } 

    U[i,j] <- gb*B[i,j]*W[i,j]/(W[i,j] + kb)
    I[i,j] <- gamma*H[i,j]*(B[i,j] + k*W0)/(B[i,j] + k)

    if (i < Nt){ # PREVENTING THE CALCULATION FROM OCCURING
      B[i+1,j] <- floor(B[i,j] + dt*((cb*U[i,j]) - (db*B[i,j]) + 
                                       (Db/(dx^2))*(Bleft - 2*B[i,j] + Bright)))

      W[i+1,j] <- floor(W[i,j] + dt*((I[i,j]  - U[i,j] - rw * 
                                        W[i,j]) + (Dw/(dx^2))*(Wleft - 2*W[i,j] + Wright)))

      H[i+1,j] <- floor(H[i,j] + dt*((R - I[i,j]) + (Dh/(dx^2))* 
                                       (Hleft - 2*H[i,j] + Hright)))
    }
  }
}