###########################################################
#                                                         #                
#                        Ú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
# (okrajové úlohy), rootSolve, PBSddesolve, sde, pomp a další.

install.packages("deSolve")
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áno 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 neomezené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:100)/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")
}

