######################################
#                                    #
#  Model dravec-kořist Gauseho typu  #
#                                    #
######################################

rk<-2
Kk<-1000
S<-2
sigma<-400
delta<--3.5   # delta = rp
c<-3

p<-function(Nk,S,sigma) {
  S*Nk/(Nk+sigma)
}

plot(c(min(casy),max(casy)),c(0,1400),cex=0,main="Model dravec-kořist Gauseho typu",xlab="Čas",ylab="Velikost populace")

for(i in 1:3) {
  N0k<-500*i
  N0p<-500*(4-i)
  
  podminky<-c(Nk=N0k,Np=N0p)
  casy<-c(0:200)/10
  soustava<-function(t,prom,param) {
    return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-p(prom[1],param[5],param[6])*prom[2],
                  prom[2]*(param[3]+param[4]*p(prom[1],param[5],param[6])))))
  }
  parametry<-c(rk,Kk,delta,c,S,sigma)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,"time"],reseni[,"Nk"],lwd=3,lty=1)
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2)  

}

## Vyšetření řešení modelu
# Np*(delta+c*p(Nk,S,sigma)) = 0
# 1. Np = 0 (predátor nikdy nebyl přítomen)
# 2. delta = - c*p(Nk, S, sigma)

library(utils)

# Ustálená velikost populace kořisti
fce<-function(Nk,delta,c,S,sigma) {
  -delta-c*S*Nk/(Nk+sigma)
}
Nklim<-uniroot(fce,c(0,1000),tol=0.1,delta=delta,c=c,S=S,sigma=sigma)$root

# Ustálená velikost populůace predátora
fce<-function(Np,Nklim,rk,Kk,S,sigma) {
  rk*Nklim*(1-Nklim/Kk)-Np*S*Nklim/(Nklim+sigma)
}
Nplim<-uniroot(fce,c(0,1000),tol=0.1,rk=rk,Kk=Kk,S=S,sigma=sigma,Nklim=Nklim)$root

lines(c(min(casy),max(casy)),rep(Nklim,2),col="red",lty=1)
lines(c(min(casy),max(casy)),rep(Nplim,2),col="red",lty=2)

















