206 lines
5.5 KiB
R
Executable File
206 lines
5.5 KiB
R
Executable File
# 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")
|