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

library("deSolve")

rk<-2
Kk<-1000
S<-2
sigma<-400
delta<--3.5   # delta = rp (musí být záporná)
c<-3

c*S*(Kk-sigma)/(Kk+sigma)

# predační (trofická) funkce Holling II
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,col="black")
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2,col="green")  

}



## Vyšetření řešení modelu (první varianta pomocí balíku nleqslv)
# 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("nleqslv")

# Ustálená velikost populace kořisti (z rovnice pro Np)
fce.k<-function(Nk) {
  delta+c*S*Nk/(Nk+sigma)
}
Nklim<-nleqslv(100,fce.k,control=list(ftol=0.0001))$x

# Ustálená velikost populace predátora (dosazením do rovnice pro Nk)
fce.p<-function(Np) {
  rk*Nklim*(1-Nklim/Kk)-Np*S*Nklim/(Nklim+sigma)
}
Nplim<-nleqslv(100,fce.p,control=list(ftol=0.0001))$x





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

#install.packages("utils")
library("utils")

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

# Ustálená velikost populace predátora (dosazením do rovnice pro Nk)
fce.p<-function(Np,Nklim,rk,Kk,S,sigma) {
  rk*Nklim*(1-Nklim/Kk)-Np*S*Nklim/(Nklim+sigma)
}
Nplim<-uniroot(fce.p,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)




## Varianty chování modelu
# 1) -delta > c * S -> populace dravce vymře
# 2) -delta < c * S -> populace koexistují
# 3) -delta > c * S * (Kk - sigma) / (Kk + sigma) -> populace se ustálí
# 3) -delta < c * S * (Kk - sigma) / (Kk + sigma) -> populace kolísají

