library("deSolve")

a1<-3
a2<-2
a3<--1  

b11<-+0
b12<-+0
b13<--1/200
b21<-+0
b22<-+0
b23<--1/300
b31<-+1/400
b32<-+1/600
b33<-+0
  
N10<-100
N20<-200
N30<-300

############################################
# Klasické řešení jako soustava tří rovnic #
############################################
podminky<-c(N1=N10,N2=N20,N3=N30)
casy<-c(0:1000)/10
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1]+param[ 4]*prom[1]+param[ 5]*prom[2]+param[ 6]*prom[3]),
                prom[2]*(param[2]+param[ 7]*prom[1]+param[ 8]*prom[2]+param[ 9]*prom[3]),
                prom[3]*(param[3]+param[10]*prom[1]+param[11]*prom[2]+param[12]*prom[3]))))
}
parametry<-c(a1,a2,a3,b11,b12,b13,b21,b22,b23,b31,b32,b33)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1200),xlab="Čas",ylab="Počet jedinců",main="Lotkův-Volterrův systém pro tři populace")
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)


#####################
# Maticová soustava #
#####################
A<-c(a1,a2,a3)
B<-matrix(c(b11,b21,b31,b12,b22,b32,b13,b23,b33),nrow=3,ncol=3)
N0<-c(N10,N20,N30)

podminky<-c("N"=N0)
casy<-c(0:1000)/10
soustava<-function(t,prom,param) {
  return(list(diag(prom)%*%(param[["A"]]+param[["B"]]%*%prom)))
}
parametry<-list("A"=A,"B"=B)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1200),xlab="Čas",ylab="Počet jedinců",main="Lotkův-Volterrův systém pro tři populace")
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)


#######################################################
# Maticová soustava alternativně (s názvy proměnných) #
#######################################################
A<-c(a1,a2,a3)
B<-matrix(c(b11,b21,b31,b12,b22,b32,b13,b23,b33),nrow=3,ncol=3)
N0<-c("N1"=N10,"N2"=N20,"N3"=N30) # podmínky musí být povinně vektor (požadavek funkce ode), nemohou být seznam

podminky<-N0
casy<-c(0:1000)/10
soustava<-function(t,prom,param) {
  return(list(diag(as.numeric(prom[names(N0)]))%*%(param[["A"]]+param[["B"]]%*%prom[names(N0)]))) #pojmenování je možné názvy prvků ve vektrou N0 (names(N0))
}
parametry<-list("A"=A,"B"=B)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1200),xlab="Čas",ylab="Počet jedinců",main="Lotkův-Volterrův systém pro tři populace")
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)


######################################
# Splněná domácí úloha (podzim 2023) #
######################################
a1<-3
a2<-3
a3<-3  

b11<--1/200
b12<-+0
b13<--1/100
b21<--1/100
b22<--1/200
b23<-+0
b31<-+0
b32<--1/100
b33<--1/200

N10<-100
N20<-200
N30<-300

A<-c(a1,a2,a3)
B<-matrix(c(b11,b21,b31,b12,b22,b32,b13,b23,b33),nrow=3,ncol=3)
N0<-c("N1"=N10,"N2"=N20,"N3"=N30) # podmínky musí být povinně vektor (požadavek funkce ode), nemohou být seznam

podminky<-N0
casy<-c(0:400)/10
soustava<-function(t,prom,param) {
  return(list(diag(as.numeric(prom[names(N0)]))%*%(param[["A"]]+param[["B"]]%*%prom[names(N0)]))) #pojmenování je možné názvy prvků ve vektrou N0 (names(N0))
}
parametry<-list("A"=A,"B"=B)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,400),xlab="Čas",ylab="Počet jedinců",main="Lotkův-Volterrův systém pro tři populace")
lines(reseni[,1],reseni[,3],col="green",lty=1,lwd=3)
lines(reseni[,1],reseni[,4],col="brown",lty=1,lwd=3)