ModeloSirdItaly-17-02-2021, R

ModeloSirdItaly-17-02-2021

# Ajuste de curvas por mínimos quadrados.

# Dados os vetores x e y, supomos que y_i=f(x_i)+erro_i, como erro_i um erro "pequeno".

# 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); n

# Suponha que S'=-aSI, I'=aSI-bI-cI; R'=bI; D'=cI. Descobrir aproximações para a, b e c.
# 
# Aproximação S_(i+1)-S_i=hS';.... Tome y=(S,I,R,D), f(y)=y-(aSI,bI+cI-aSI,bI,cI)
# 
# Para determinar os valores dos parâmetros resolvemos o problema de minimização 
#
#        min g(u)
# 
# em que g(u)= sum (y_(i+1)-f(x_i))^2, sendo u=(a,b,c).
# 
# Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial 
#
#      u'(t)=-grad g(u(t))
#      u(0)=u_0

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

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

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

#---- Início da rotin para o método de Euler.

ZeroEuler<-function(u0,t,m){
w=u0
h=t/m
for (i in 1:m){w=w-h*gradng(w)}
w
}

u0=c(1,1,1)

print("Teste condição inicial"); g(u0)
u=ZeroEuler(u0,20,20000)
#---------------------------------------------------------------------------------------

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
fap<-function(s){a/(b+c*exp(-a*(s-1930)/10))}  # Aproximação da função.

print("Gráfico para ajuste")
curve(fap,1930,2010)                           # Teste de ajuste
points(seq(1930,2010,by=10),y,col="red")

print("População estimada em 2021");fap(2021)

# plot(tt,mf(Y[1,],Y[2,]),col = "blue",'l') # curva f(u(t),v(t)) para caminho de maximização obtido pelo método de Euler.


run \| edit \| history \| help	0