###########################################################
#                                                         #                
#                        Úvod do R                        #
#        Bi3101 Úvod do matematického modelování          #
#                                                         #
#  V následujícím skriptu byly využity úryvky z článku:   #
#           Solving Differential Equations in R           #
# (Karline Soetaert, Thomas Petzoldt a R. Woodrow Setzer) #
#                                                         #
###########################################################

# Existuje více různých balíků pro numerické řešení (soustav) diferenciálních rovnic v R.
# Pro řešení obyčejných diferenciálních rovnic (tj. takových, kde je jedna nezávisla
# proměnná) lze výhodně použít balík deSolve. Další jsou např. balíky ReacTran pro
# parciální differenciální rovnice, bvpSolve pro soustavy s okrajovými podmínkami
# (okrajvé úlohy), rootSolve, PBSddesolve, sde, pomp a další.

library("deSolve")

# Základní funkcí pro (numerické) řešení (soustav) diferenciálních rovnic je funkce
# ode, jejíž syntaxe je ode(podminky,soustava,casy,parametry,metoda):
#
# podminky………počáteční podmínky ve formě vektoru rovnic
# casy…………………časy, ve kterých bude získánio numerické řešení, první hodnota je počátek
# soustava………funkce vracející seznam derivací jednotlivých proměnných
# parametry……další případné parametry funkce soustava
# metoda……………metoda použitá pro řešení, implicitně "lsoda", která umí vybrat jinou metodu

# Pro zápis rovnice neomzeného růstu populace lze psát:
N0<-100
a <-7/12
b <-1/(11*12)


# Diskrétní případ

podminky<-c(N=N0)
casy<-0:10
soustava<-function(t,prom,param) {
  return(list(c(prom[1]+prom[1]*(param[1]-param[2]))))
}
parametry<-c(a,b)
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:10
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]-param[2]))))
}
parametry<-c(a,b)
metoda<-"lsoda"

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


#######################
#                     #
#  Modifikace modelu  #
#                     #
#######################

# Stejně jako v předchozí práci v Maple, i nyní zavedem kvadratický člen,
# který zbrzdí neomezený rozvoj populace - úživnost prostředí.
N0<-20
a <-7/12
b <-1/(11*12)
c <-1/(11*12)


# Diskrétní případ

podminky<-c(N=N0)
casy<-0:12
soustava<-function(t,prom,param) {
  return(list(c(prom[1]+prom[1]*(param[1]-param[2])-prom[1]^2*param[3])))
}
parametry<-c(a,b,c)
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:12
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]-param[2])-prom[1]^2*param[3])))
}
parametry<-c(a,b,c)
metoda<-"lsoda"

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


# Vykreslení pro 15 různých počátečních podmínek

casy<-0:10
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]-param[2])-prom[1]^2*param[3])))
}
parametry<-c(a,b,c)
metoda<-"lsoda"

plot(c(0,10),c(0,150),cex=0,xlab="t",ylab="N")
for (i in 1:15) {
  podminky<-c(N=10*i)
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,1],reseni[,2],lwd=3,col="brown")
}


# Vykreslení pro 15 různých počátečních podmínek, hladší křivky

casy<-seq(0,10,0.1)
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]-param[2])-prom[1]^2*param[3])))
}
parametry<-c(a,b,c)
metoda<-"lsoda"

plot(c(0,10),c(0,150),cex=0,xlab="t",ylab="N")
for (i in 1:15) {
  podminky<-c(N=10*i)
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,1],reseni[,2],lwd=3,col="brown")
}


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

# 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)
  }
}


# 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:12
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")

###########################
#                         #
#  Interagující populace  #
#                         #
###########################


# Model konkurekce

r1<-1
r2<-0.6
K1<-1000
K2<-1500
C1<-600
C2<-1000
gamma1<-0.5
gamma2<-0.3
N10<-300
N20<-300

kappa<-function(N,K,C,gamma){
  return((K+C*gamma*N)/(1+gamma*N))
}

podminky<-c(N1=N10, N2=N20)
casy<-0:15
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1])*(1-(prom[1]/kappa(prom[2],param[3],param[5],param[7]))),
                prom[2]*(param[2])*(1-(prom[2]/kappa(prom[1],param[4],param[6],param[8]))))))
}
parametry<-c(r1,r2,K1,K2,C1,C2,gamma1,gamma2)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1100))
lines(reseni[,1],reseni[,3],col="green",lty=2,lwd=3)

# vnořené cykly
plot(c(0,15),c(300,1200),cex=0,xlab="t",ylab="N")
for(i in 1:4){
  for(j in 1:4){
    N10<-300*i
    N20<-300*j
    podminky<-c(N1=N10, N2=N20)
    reseni<-ode(podminky,casy,soustava,parametry,metoda)
    lines(reseni[,1],reseni[,2],lwd=3,col="black")
    lines(reseni[,1],reseni[,3],lwd=3,col="green",lty=2)
  }
}


# Model symbiózy

r1<-1
r2<-0.6
K1<-1000
K2<-1500
C1<-1600
C2<-2000
gamma1<-0.5
gamma2<-0.3
N10<-300
N20<-300

kappa<-function(N,K,C,gamma){
  return((K+C*gamma*N)/(1+gamma*N))
}

podminky<-c(N1=N10, N2=N20)
casy<-0:15
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1])*(1-(prom[1]/kappa(prom[2],param[3],param[5],param[7]))),
                prom[2]*(param[2])*(1-(prom[2]/kappa(prom[1],param[4],param[6],param[8]))))))
}
parametry<-c(r1,r2,K1,K2,C1,C2,gamma1,gamma2)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1100))
lines(reseni[,1],reseni[,3],col="green",lty=2,lwd=3)

# vnořené cykly
plot(c(0,15),c(300,2100),cex=0,xlab="t",ylab="N")
for(i in 1:4){
  for(j in 1:4){
    N10<-500*i
    N20<-600*j
    podminky<-c(N1=N10, N2=N20)
    reseni<-ode(podminky,casy,soustava,parametry,metoda)
    lines(reseni[,1],reseni[,2],lwd=3,col="black")
    lines(reseni[,1],reseni[,3],lwd=3,col="green",lty=2)
  }
}