########################
#                      #
#  Model soudného dne  #
#                      #
########################

# Model exponenciálního růstu dN(t)/d(t) = r ⋅ N(t) vedl v našich předchozích úvahách k nerealistickému, extrémnímu růstu.
# Při pokusech o modelování reálné světové lidské populace nicméně i tento model vykazuje stále příliš pomalý růst.
# Modeláři zjistili, že vhodnější je použití ještě strmějšího modelu, zvaného „model soudného dne“: dN(t)/d(t) = r ⋅ N(t)^(1+c).
# Vyřešte obecně, na jaké hodnotě velikosti populace se daný model ustálí a okomentujte, proč.

library("deSolve")
library("nleqslv")

r<-0.5
c<-0.01
N0<-100


podminky<-c(N=N0)
casy<-(0:100)/10
soustava<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]^(1+param[2]))))
}
parametry<-c(r,c)
metoda<-"lsoda"

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

# Kdy nastane soudný den

#dN(t)/d(t) = r ⋅ N(t)^(1+c)       #řešíme diferenciální rovnici

#dN(t)/ N(t)^(1+c) = r ⋅ d(t)      #převedení diferenciálů na protilehlé strany rovnice

#dN(t) ⋅ N(t)^(-1-c) = r ⋅ d(t)    #záporný exponent ze jmenovatele do čitatele

#N(t)^(-c) / (-c) = K + r ⋅ t          #zintegrujeme obě strany rovnice

#c / N(t)^(-c) = -1 / (K + r ⋅ t)

#N(t)^(c) = -1 / (c ⋅ (K + rt))

#N(t) = (-1 / (c ⋅ (K + rt)))^(1/c)

# Dále spočítáme vztah konstanty K a počátečního počtu N0.

#N(t)^(-c) / (-c) = K + r ⋅ t 

#N0^(-c) / (-c) = K + r ⋅ 0

#N0^(-c) / (-c) = K

# A dosadíme

#N(t) = (-1 / (c ⋅ (N0^(-c) / (-c) + rt)))^(1/c)

#N(t) = (-1 / ((-N0^(-c) + crt)))^(1/c)

#N(t) = (1 / ((N0^(-c) - crt)))^(1/c)

#N(t) = (N0^c / (1 - crt ⋅ N0^c))^(1/c)

#N(t) = N0 / ((1 - crt ⋅ N0^c)^(1/c))

# Obecně tedy platí, že N(t) bude nekonečné, když bude jmenovatel na pravé straně nulový

# (1 - crt ⋅ N0^c)^(1/c) = 0

# 1 - crt ⋅ N0^c = 0

# 1 = crt ⋅ N0^c

# t = 1 / (cr ⋅ N0^c)

# V našem případě konkrétně

1/(c*r*N0^c)

# Zaspáno jako funkce
N<-function(t) {
  N0/((1-c*r*t*N0^c)^(1/c))
}

# Vykresleno
plot(1:1000,N(1:1000),type="l",ylim=c(0,1e20),lwd=3,col="brown") # hodnoty jsou extrémně velké, pomůžeme si logaritmem
plot(1:1000,log10(N(1:1000)),type="l",ylim=c(-100,400),lwd=3,col="brown") # graficky - řešení je kolem hodnoty 190
plot(1:1000,-log10(log10(N(1:1000))),type="l",ylim=c(-2.5,1),lwd=3,col="brown") # dvakrát zlogaritmováno a obráceno, pro přesnou identifikaci „špičky“, takový úlet

# Zkusíme vyřešit též numericky - protože příkaz nleqslv() umí najít kořen (nulové řešení) a my hledáme nekonečné, obrátíme hodnotu N
obracene.N<-function(t) {
  1/N(t)
}

soudnyden.start<-100
nleqslv(soudnyden.start,obracene.N)[[1]]                                     #změny jsou tak malé, že bez nastavení citlivosti se vůbec nehne z iniciální hodnoty
nleqslv(soudnyden.start,obracene.N,control=list(ftol=1e-161,maxit=1e6))[[1]] #citlivost ftol zvýšená na maximální povolenou hodnotu 1e-161

# Přesnost můžeme ještě zvýšit vynásobením nějakým velikým číslem
obracene.N.zvetsene<-function(t) {
  1e60/N(t)
}

soudnyden.start<-100
nleqslv(soudnyden.start,obracene.N.zvetsene,control=list(ftol=1e-161,maxit=1e6))[[1]] 







































