# Modelo SIR
library(deSolve)
#tamaño poblacional
N = 1
#estado inicial de los compartimentos
init <- c(S = 1-1e-6,
I = 1e-6,
R = 0)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 2.0,
gamma = 0.2)
#crear la función con las ODE
sir <- function(times, init, param) {
with(as.list(c(init, param)), {
#ecuaciones diferenciales
dS <- -beta * S * I
dI <-  beta * S * I - gamma * I
dR <-                 gamma * I
#resultados de las tasas de cambio
return(list(c(dS, dI, dR)))
})
}
#intervalo de tiempo y resolución
times <- seq(0, 20, by = 1)
#resolver el sistema de ecuaciones con función 'ode'
out <- ode(y = init, times = times, func = sir, parms = param)
#cambiar out a un data.frame
out <- as.data.frame(out*N) #aqui puede multiplicar 'out' por N
# ver resultados
head(out)
plot(out)
plot(out$time,out$S)
line(out$I)
abline(out$I)
lines(out$I)
lines(out$R)
plot(out$time,out$S,type="b")
lines(out$I)
lines(out$R)
plot(out$time,out$S,type="l")
lines(out$I)
lines(out$R)
plot(out$time,out$S,type="l")
lines(out$I,color="red")
lines(out$R,color="blue")
plot(out$time,out$S,type="l")
lines(out$I,col="red")
lines(out$R,col="blue")
plot(out$time,out$S,type="l",labx="Tiempo",laby="Población")
lines(out$I,col="red")
lines(out$R,col="blue")
plot(out$time,out$S,type="l",xlab="Tiempo",ylab="Población")
lines(out$I,col="red")
lines(out$R,col="blue")
Data <- read.delim("D:/Cursos/EDO/Data.txt")
View(Data)
#exponecial
r=0.15
Y<-seq(1:50)
t<-seq(1:50)
for (i in 1:49) {Y[i+1]<r*exp(Y[i])}
plot(t,Y)
#exponecial
r=1.15
Y<-seq(1:50)
t<-seq(1:50)
for (i in 1:49) {Y[i+1]<r*exp(Y[i])}
plot(t,Y)
#exponecial
r=1.15
Y<-seq(1:50)
t<-seq(1:50)
for (i in 1:49) {Y[i+1]<r*exp(Y[i])}
plot(t,Y)
# Logistico
Dias<-c(105,112,119,126,133,140,147,154,161,168,175)
Peso<-c(59.4,66.1,72.8,79.6,86.4,93.2,99.8,106.4,112.8,119.2,125.3)
#Regresión no lineal
Reg<-nls(Peso~k/(1+b*exp(-r*Dias)),star=c(k=80,b=0.1, r=1),control=nls.control(warnOnly=TRUE))
summary(Reg)
K=180
b1=513
r=0.72
Y<-seq(1:50)
t<-seq(1:50)
G<-seq(1:50)
for (i in 1:49) {Y[i+1]<-K/(1+b1*exp(-r*i))}
plot(t,Y,xlab="Tiempo (semanas)",ylab="Peso Vivo")
for (i in 1:49) {G[i+1]<-(Y[i+1]-Y[i])}
plot(t,G)
# Presa-depredador Euler
a<-0.03
b<-0.0009
c<-0.0005
d<-0.05
N<-1000
x<-seq(1:N)
y<-seq(1:N)
paso<-1.5
x[1]<-100
y[1]<-100
for (i in 1:N) {x[i+1]<-x[i]+paso*(a*x[i]-b*x[i]*y[i])
y[i+1]<-y[i]+paso*(c*x[i]*y[i]-d*y[i])}
plot(x,ylim=c(0,max(x+y)),type='l',ylab="Densidad Poblacional")
lines(y,lwd=2,col="red")
# Modelo SIR
library(deSolve)
#tamaño poblacional
N = 1
#estado inicial de los compartimentos
init <- c(S = 1-1e-6,
I = 1e-6,
R = 0)
#parámetros del modelo (coeficientes de las variables)
param <- c(beta = 2.0,
gamma = 0.2)
#crear la función con las ODE
sir <- function(times, init, param) {
with(as.list(c(init, param)), {
#ecuaciones diferenciales
dS <- -beta * S * I
dI <-  beta * S * I - gamma * I
dR <-                 gamma * I
#resultados de las tasas de cambio
return(list(c(dS, dI, dR)))
})
}
#intervalo de tiempo y resolución
times <- seq(0, 20, by = 1)
#resolver el sistema de ecuaciones con función 'ode'
out <- ode(y = init, times = times, func = sir, parms = param)
#cambiar out a un data.frame
out <- as.data.frame(out*N) #aqui puede multiplicar 'out' por N
# ver resultados
head(out)
plot(out$time,out$S,type="l",xlab="Tiempo",ylab="Población")
lines(out$I,col="red")
lines(out$R,col="blue")
N<-1000
t0<-0
tN<-600
h<-(tN-t0)/N
alpha1<-1/14
alpha2<-1/6
alpha3<-1/4
f<-30
miu1<-0.0025
miu2<-0.004
miu3<-0.0028
ro<-0.0068
pi<-50
L<-0.00165
d<-0.008
k<-0.0008
H<-c(1:N)
J<-c(1:N)
A<-c(1:N)
R<-c(1:N)
for (i in 1:N) {H[i+1]<-H[i]+h*(alpha3*f*A[i]*R[i]-(alpha1+miu1)*J[i])
J[i+1]<-J[i]+h*(alpha1*H[i]-(alpha2+miu2)*J[i])
A[i+1]<-A[i]+h*(alpha2*J[i]*R[i]-alpha3*d*A[i])
R[i+1]<-R[i]+h*(ro*R[i]*(1-pi*R[i])-L*J[i]*R[i]-d*A[i]*R[i]-k*R[i])}
N<-H+J+A
plot(N,ylim=c(0,max(N)),type='l',xlab="Tiempo(horas)",ylab="Densidad Poblacional")
lines(H,lwd=2,col="green")
lines(J,lwd=2,col="brown")
lines(A,lwd=2,col="blue")
lines(R,lwd=2,col="orange")
title("Dinámica poblacional del nematodo en la raíz")
text(1020,0.05, labels = "R", col = "orange")
text(1020,13, labels = "A", col = "blue")
text(1020,5, labels = "J", col = "brown")
text(1020,8, labels = "H", col = "green")
text(1020,20, labels = "N", col = "black")
plot(N,ylim=c(0,max(N)),type='l',xlab="Tiempo(horas)",ylab="Densidad Poblacional")
lines(H,lwd=2,col="green")
lines(J,lwd=2,col="brown")
lines(A,lwd=2,col="blue")
lines(R,lwd=2,col="orange")
title("Dinámica poblacional del nematodo en la raíz")
text(1020,20, labels = "R", col = "orange")
text(1020,25, labels = "A", col = "blue")
text(1020,40, labels = "J", col = "brown")
text(1020,50, labels = "H", col = "green")
text(1020,60, labels = "N", col = "black")
plot(N,ylim=c(0,max(N)),type='l',xlab="Tiempo(horas)",ylab="Densidad Poblacional")
lines(H,lwd=2,col="green")
lines(J,lwd=2,col="brown")
lines(A,lwd=2,col="blue")
lines(R,lwd=2,col="orange")
title("Dinámica poblacional del nematodo en la raíz")
text(800,20, labels = "Raiz", col = "orange")
text(800,30, labels = "Adulto", col = "blue")
text(800,40, labels = "Juveniles", col = "brown")
text(800,50, labels = "Huevos", col = "green")
text(800,60, labels = "Total", col = "black")
plot(N,ylim=c(0,max(N)),type='l',xlab="Tiempo(horas)",ylab="Densidad Poblacional")
lines(H,lwd=2,col="green")
lines(J,lwd=2,col="brown")
lines(A,lwd=2,col="blue")
title("Dinámica poblacional del nematodo en la raíz")
text(800,30, labels = "Adulto", col = "blue")
text(800,40, labels = "Juveniles", col = "brown")
text(800,50, labels = "Huevos", col = "green")
text(800,60, labels = "Total", col = "black")