SimulacaoItaly-17-02-2021-linearizado, R

Language: Layout:

##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

# Seguem os comando para as definições de funções necessárias.

g<-function(u){    # função g(u)
a=u[1];b=u[2]; c=u[3]
y=matrix(0,4,n)
y[1,]=S; y[2,]=I; y[3,]=R; y[4,]=D
p=c(0,0,0,0)
for (i in 1:(n-1)){
p=p+( y[,i+1]-(y[,i]-c(a*S[i]*I[i],b*I[i]+c*I[i]-a*S[i]*I[i],b*I[i],c*I[i])) )^2
}
sum(p)/2
}

# Teste do código para g(u).
print(" Teste g"); g(c(1,2,1))

v=matrix(0,4*(n-1),1)
A=matrix(0,4*(n-1),3)
for ( i in 1:(n-1)){
B=matrix(0,4,3); B[1,]=c(-S[i]*I[i],0,0); B[2,]=c(S[i]*I[i],I[i],I[i]); B[3,]=c(0,I[i],0); B[4,]=c(0,0,I[i])
A[(4*i-3):(4*i),]=B
v[(4*i-3):(4*i)]=c(S[i+1]-S[i],I[i+1]-I[i],R[i+1]-R[i],D[i+1]-D[i])
}

gradng<-function(u){   # função gradiente de g(u)=|Au-v|^2.
a=u[1];b=u[2]; c=u[3]
t(A)%*%(A%*%u-v)
}

# Teste da função gradg(u)
print(" Teste gardiente de g"); gradng(c(1,2,1))

u=solve(t(A)%*%A,t(A)%*%v)
print("Soulução"); u                                
a=u[1]; b=u[2]; c=u[3]

print(" Resíduo da aproximação"); g(u)
print("Gradiente no ponto aproximado"); gradng(u)
print("Parâmetros a, b e c encontrados"); a;b;c
# Teste ajuste

RungK4v<-function(f,u0,n){m=length(u0); Y=matrix(0,m,n+1); Y[,1]=u0; 
for (i in 1:n){ 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/100; h
times = seq(0, b, by = h); n=length(times)-1; n

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')

[ + ] Show input
Absolute running time: 1.33 sec, cpu time: 1.58 sec, memory peak: 55 Mb, absolute service time: 1,57 sec