######################################################
#                                                    #
#  Populace pod tlakem nespecializovaného predátora  #
#                                                    #
######################################################

library("deSolve")

# Definice predační funkce p (Nk = kritické N, nad kterým je hodnota p rovna S = saturace)

p<-function(N,Nk,S) {
  if (N <= Nk) {
    return((S/Nk)*N)
  }
  else{
    return(S)
  }
}

q<-function(N,Nk,S) {
  return((S*N)/(N+Nk))
}

r<-function(N,Nk,S) {
  return(2*S/(1+exp(-N/1000))-S)
}

plot(c(0:10000),r(c(0:10000),Nk,S),type="l")


# Diskrétní případ

N0<-500
r<-1
K<-1000
Nk<-200
S<-300

podminky<-c(N=N0)
casy<-0:12
soustava<-function(t,prom,param) {
  return(list(c(prom[1]+param[1]*prom[1]*(1-prom[1]/param[2])-p(prom[1],param[3],param[4]))))
}
parametry<-c(r,K,Nk,S)
metoda<-"iteration"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="p",col="brown",lwd=3,which="N")


# Spojitý případ

podminky<-c(N=N0)
casy<-0:100
soustava<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-p(prom[1],param[3],param[4]))))
}
parametry<-c(r,K,Nk,S)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="brown",lwd=3,which="N",ylim=c(0,1000))

soustava2<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-q(prom[1],param[3],param[4]))))
}

reseni2<-ode(podminky,casy,soustava2,parametry,metoda)
lines(reseni2[,1],reseni2[,2],col="cornflowerblue",lwd=3)

soustava3<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-r(prom[1],param[3],param[4]))))
}

reseni3<-ode(podminky,casy,soustava3,parametry,metoda)
lines(reseni3[,1],reseni3[,2],col="gold",lwd=3)


# Analýza řešení

r*N*(1-N/K)-p(N)=0

a^2+3*a-0.5*b=b

leftside<-function(a) a^2+3*a-0.5*b
b<-leftside(4)
