Использование событий в deSolve для предотвращения отрицательных переменных состояния, R

Question

Я моделирую изменение численности в пищевой сети видов, используя ODE и deSolve в R. Очевидно, что популяции не должны быть меньше нуля. поэтому я добавил функцию события и запустил ее, как показано ниже. хотя ответы меняются с того момента, когда я использовал функцию события, но она все же производит отрицательные значения. Что не так?

#using events in a function to distinguish and address the negative abundances 
eventfun <- function(t, y, parms){ 
    y[which(y<0)] <- 0 
    return(y) 
} 

# =============================== main code 
max.time = 100 
start.time = 50 

initials <- c(N, R) 

#parms <- list(webs=webs, a=a, b=b, h=h, m=m, basals=basals, mu=mu, Y=Y, K=K, no.species=no.species, flow=flow,S=S, neighs=neighs$neighs.per, dispers.maps=dispers.maps) 
temp.abund <- ode(y=initials, func=solve.model, times=0:max.time, parms=parms, events = list(func = eventfun, time = 0:max.time)) 

и здесь функция ОДА (если это помогает в поиске проблемы):

solve.model <- function(t, y, parms){ 

y <- ifelse(y<1e-6, 0, y) 
with(parms,{ 
    # return from vector form into matrix form for calculations 
    (R <- as.matrix(y[(max(no.species)*length(no.species)+1):length(y)])) 
    (N <- matrix(y[1:(max(no.species)*length(no.species))], ncol=length(no.species))) 

    dy1 <- matrix(nrow=max(no.species), ncol=length(no.species)) 
    dy2 <- matrix(nrow=length(no.species), ncol=1) 

    no.webs <- length(no.species) 
    for (i in 1:no.webs){ 

     species <- no.species[i] 
     (abundance <- N[1:species,i]) 

     adj <- as.matrix(webs[[i]]) 
     a.temp <- a[1:species, 1:species]*adj 
     b.temp <- b[1:species, 1:species]*adj 
     h.temp <- h[1:species, 1:species]*adj 

     (sum.over.preys <- abundance%*%(a.temp*h.temp)) 
     (sum.over.predators <- (a.temp*h.temp)%*%abundance) 

     #Calculating growth of basal 
     (basal.growth <- basals[,i]*N[,i]*(mu*R[i]/(K+R[i])-m)) 

     # Calculating growth for non-basal species D 
     no.basal <- rep(1,len=species)-basals[1:species] 

     predator.growth<- rep(0, max(no.species)) 
     (predator.growth[1:species] <- ((abundance%*%(a.temp*b.temp))/(1+sum.over.preys)-m*no.basal)*abundance) 


     predation <- rep(0, max(no.species)) 
     (predation[1:species] <- (((a.temp*b.temp)%*%abundance)/t(1+sum.over.preys))*abundance) 

     (pop <- basal.growth + predator.growth - predation) 


     dy1[,i] <- pop 

     dy2[i] <- 0.0005 #to consider a nearly constant value for the resource 

    } 

    #Calculating dispersals .they can be easily replaced 
    # by adjacency maps of connections between food webs arbitrarily! 
    disp.left <- dy1*d*dispers.maps$left.immig 
    disp.left <- disp.left[,neighs[,2]] 

    disp.right <- dy1*d*dispers.maps$right.immig 
    disp.right <- disp.right[,neighs[,3]] 

    emig <- dy1*d*dispers.maps$emigration 
    mortality <- m*dy1 

    dy1 <- dy1+disp.left+disp.right-emig 

    return(list(c(dy1, dy2)))  
}) 
} 

спасибо вам за вашу помощь

Answer 1

Я имел успех, используя аналогичную функцию события определяется следующим образом:

eventfun <- function(t, y, parms){ 
    with(as.list(y), { 
    y[y < 1e-6] <- 0 
    return(y) 
    }) 
}