SimulacaoItaly-17-02-2021, R

SimulacaoItaly-17-02-2021

##Modelo covid-19 SIRD

# Resolução numérica por Runge-Kutta 4

RungK4<-function(f,u0,n){m=length(u0); h=1/n; Y=u0; 
for (i in 1:n){ k1=f(Y)
k2=f(Y+h*k1/2); k3=f(Y+h*k2/2); k4=f(Y+h*k3)
Y=Y+( k1+ 2*k2 + 2*k3+ k4)*h/6}; Y}
# Fim Runge-Kutta

sirmod <- function(y,parms) {
S = y[1]
I = y[2]
R = y[3]
D = y[4]
beta = parms[1] 
gamma = parms[2] 
mu = parms[3] 
dS = - beta * S * I
dR = gamma * I
dD = mu* I
dI = -dS-dR-dD
res = c(dS, dI, dR, dD);res
}
# Fim do modelo SIRD

# Dados Italia 
ItalyI=c(2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,20,62,155,229,322,453,655,888,1128,1694,2036,2502,3089,3858,4636,5883,7375,9172,10149,12462,12462,17660,21157,24747,27980,31506,35713,41035,47021,53578,59138,63927,69176,74386,80589,86498,92472,97689,101739,105792,110574,115242,119827,124632,128948,132547,135586,139422,143626,147577,152271,156363,159516,162488,165155,168941,172434,175925,178972,181228,183957,187327,189973)
ItalyR=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,2,1,1,3,45,46,46,83,149,160,276,414,523,589,622,724,724,1045,1045,1439,1966,2335,2749,2941,4025,4440,4440,6072,7024,7024,8326,9362,10361,10950,12384,13030,14620,15729,16847,18278,19758,20996,21815,22837,24392,26491,28470,30455,32534,34211,35435,37130,38092,40164,42727,44927,47055,48877,51600,54543,57576)
ItalyD=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,2,3,7,10,12,17,21,29,34,52,79,107,148,197,233,366,463,631,827,827,1266,1441,1809,2158,2503,2978,3405,4032,4825,5476,6077,6820,7503,8215,9134,10023,10779,11591,12428,13155,13915,14681,15362,15887,16523,17127,17669,18279,18849,19468,19899,20465,21067,21645,22170,22745,23227,23660,24114,24648,25085,25549)
pop=max(ItalyI); pop
ItalyS=pop-ItalyI
ItalyI=ItalyI-ItalyD-ItalyR
plot(1:length(ItalyS),ItalyS/pop);               S=ItalyS/pop
points(1:length(ItalyI),ItalyI/pop,col="blue");  I=ItalyI/pop
points(1:length(ItalyD),ItalyD/pop,col="red");   D=ItalyD/pop
points(1:length(ItalyR),ItalyR/pop,col="green"); R=ItalyR/pop

n=length(ItalyI)

# Inicio do ajuste

p=1
n=length(S); S=S[p:n]; I=I[p:n]; R=R[p:n]; D=D[p:n]; n=length(S)
b=length(S) ;b
h=1; h
times = seq(0, length(S) , by = h); n=length(times)-1; n

# Funcional custo
N=pop; S0= N-2 ;I0=2; R0= 0 ; D0=  0; C0=R0+I0+D0

start = c(S0,I0,R0,D0)/N

E<-function(v){
n=length(I)
parms=v
u=start
f<-function(y){sirmod(y,parms)}
p=(u[1]-S[1])^2+(u[2]-I[1] )^2 +(u[3]-R[1] )^2+(u[4] -D[1])^2
for ( j in 2:n){
u=RungK4(f,u,10)
p=p+( u[1]-S[j])^2+( u[2]-I[j] )^2 +( u[3]-R[j] )^2+( u[4] -D[j])^2
                   }
p}

hn=10^(-5)
GradE<-function(s){
u=s;v=s; w=s
u[1]=s[1]+hn; v[2]=s[2]+hn; w[3]=s[3]+hn
p=c(E(u)-E(s),E(v)-E(s),E(w)-E(s))/hn
p}

# teste

u0=c(1,1,1)*0

E(u0); GradE(u0)

m=10^4

he=1

#-- Inicio iteração
u=u0

for ( i in 1:m){ 
v=u-he*GradE(u)
if (E(v)<E(u)){u=v;he=10*he}
else {he=he/10}
}
u

E(u)

GradE(u)

# Teste ajuste 
RungK4v<-function(f,u0,nn){m=length(u0); Y=matrix(0,m,nn+1); Y[,1]=u0; 
for (i in 1:nn){ t=i*h; k1=f(Y[,i])
k2=f(Y[,i]+h*k1/2); k3=f(Y[,i]+h*k2/2); k4=f(Y[,i]+h*k3)
Y[,i+1]=Y[,i]+( k1+ 2*k2 + 2*k3+ k4)*h/6}; Y}

parms=u

# Dados para simular 
b=length(ItalyI) ;b
h=1/10; h
times = seq(0, b, by = h); nn=length(times)-1; nn

N=pop; S0= N-2 ;I0=2; R0= 0 ; D0=  0; C0=R0+I0+D0
start = c(S0,I0,R0,D0)/N

f<-function(y){sirmod(y,parms)}

out=RungK4v(f,start,length(times)-1)
##plotar o gráfico de S, para plotar E e I, basta trocar por y=out$E ou y=out$I
points(times, out[1,], ylab="S" , xlab="t" , type="l")
points(times, out[2,], ylab="I" , xlab="t" , type="l", col= 'blue')
points(times, out[3,], ylab="R" , xlab="t" , type="l", col= 'green')
points(times, out[4,], ylab="D" , xlab="t" , type="l", col= 'red')


run \| edit \| history \| help	0