Estudo do Impacto de um Choque Epidêmico

This document is a product of my work as student at São Carlos Federal University (UFSCar), Sorocaba, 2024.

true

Instalação dos Pacotes e Checagem do Diretório

Show code
#install.packages("pacman")
  library("pacman")
  p_load(phaseR, knitr, deSolve)

Desenvolvimento do Modelo

Show code
#Introdução de um indivíduo infectado (x) à uma população de n pessoas saudáveis (y), sujeitos a uma taxa de recuperação e uma taxa de mortalidade

#  y + x = n + 1 --> Isolando o y --> y = n + 1 - x

#E.D.O. --> y' = (y + x) - [k * y * (x)] --> y' = (n+1) - [k * y * (n + 1 - y)]

#r <- tr*((n+1-(yt-1) + tc*((yt-1))
#x <- ((n+1-(yt-1) + tc*((yt-1))
#m <- tm*((n+1-(yt-1) + tc*((yt-1))
#y <- nv + r - x
#nvt = nvt-1 - m

#tc = 0.7 #taxa de contágio
#n = 1000 #população total em t = 0
#y0 = 1000 ##população saudável em t=0
#tm = 0.2 #taxa de mortalidade
#tr = 0.8 #taxa de recuperação
#x0 = 1 #população infectada em t=0
#m0 = 0 #população morta em t=0
#nv0 <- nv0 - m0 #população viva em t=0 => nv=0

Cálculo e Demonstração Gráfica

Show code
# Parâmetros
  tc <- 0.7
  n <- 1000
  y0 <- 1000
  tm <- 0.2
  tr <- 0.8
  x0 <- 1
  m0 <- 0
  nv0 <- n - m0

# Função epidemiológica
  epidem <- function(t, y, parameters){
    tr <- parameters$tr
    nv <- parameters$nv
    x <- parameters$x
    dy <- nv+tr-x
    list(dy)
  }

# Vetores para armazenar os resultados
  x <- numeric(40)
  r <- numeric(40)
  m <- numeric(40)
  nv <- numeric(40)
  y <- numeric(40)

# Loop para as próximas 40 iterações
  for (i in 1:40) {
    if (i > 1) {
       nv[i] <- nv[i - 1] - m[i - 1]
        x[i] <- (nv[i] + 1 - y[i - 1] + tc * y[i - 1])
        r[i] <- tr * x[i - 1]
        m[i] <- tm * x[i - 1]
        y[i] <- nv[i] + r[i] - x[i]
    } else {
      # Condições iniciais para a primeira iteração
        nv[1] <- nv0 - m0
        x[1] <- (nv[1] + 1 - y0 + tc * y0)
        r[1] <- tr * x0
        m[1] <- tm * x0
        y[1] <- nv[1] + r[1] - x[1]
    }
  }

# Exibindo os resultados
  resultados <- data.frame(Periodo = 1:40, x, r, m, nv, y)
  print(resultados)
   Periodo           x           r           m           nv
1        1 701.0000000   0.8000000   0.2000000 1.000000e+03
2        2 910.8600000 560.8000000 140.2000000 9.998000e+02
3        3 665.6780000 728.6880000 182.1720000 8.596000e+02
4        4 401.6450000 532.5424000 133.1356000 6.774280e+02
5        5 302.7947800 321.3160000  80.3290000 5.442924e+02
6        6 296.1193140 242.2358240  60.5589560 4.639634e+02
7        7 281.3804710 236.8954512  59.2238628 4.034044e+02
8        8 237.5047539 225.1043768  56.2760942 3.441806e+02
9        9 189.3704258 190.0038032  47.5009508 2.879045e+02
10      10 154.8421769 151.4963406  37.8740852 2.404035e+02
11      11 132.4121411 123.8737415  30.9684354 2.025295e+02
12      12 114.3637002 105.9297129  26.4824282 1.715610e+02
13      13  97.1404790  91.4909602  22.8727400 1.450786e+02
14      14  81.3771268  77.7123832  19.4280958 1.222058e+02
15      15  68.2154205  65.1017015  16.2754254 1.027778e+02
16      16  57.6031165  54.5723364  13.6430841 8.650233e+01
17      17  48.8177783  46.0824932  11.5206233 7.285924e+01
18      18  41.3014318  39.0542226   9.7635557 6.133862e+01
19      19  34.8476403  33.0411454   8.2602864 5.157506e+01
20      20  29.3842064  27.8781122   6.9695281 4.331478e+01
21      21  24.8026440  23.5073651   5.8768413 3.634525e+01
22      22  20.9534166  19.8421152   4.9605288 3.046841e+01
23      23  17.7007469  16.7627332   4.1906833 2.550788e+01
24      24  14.9462359  14.1605975   3.5401494 2.131720e+01
25      25  12.6175789  11.9569887   2.9892472 1.777705e+01
26      26  10.6528621  10.0940631   2.5235158 1.478780e+01
27      27   8.9955831   8.5222897   2.1305724 1.226428e+01
28      28   7.5964138   7.1964665   1.7991166 1.013371e+01
29      29   6.4144650   6.0771310   1.5192828 8.334594e+00
30      30   5.4161333   5.1315720   1.2828930 6.815311e+00
31      31   4.5731933   4.3329066   1.0832267 5.532418e+00
32      32   3.8615521   3.6585546   0.9146387 4.449192e+00
33      33   3.2606947   3.0892417   0.7723104 3.534553e+00
34      34   2.7533126   2.6085558   0.6521389 2.762243e+00
35      35   2.3248579   2.2026500   0.5506625 2.110104e+00
36      36   1.9630723   1.8598863   0.4649716 1.559441e+00
37      37   1.6575930   1.5704579   0.3926145 1.094470e+00
38      38   1.3996547   1.3260744   0.3315186 7.018550e-01
39      39   1.1818540   1.1197238   0.2799309 3.703364e-01
40      40   0.9979436   0.9454832   0.2363708 9.040549e-02
              y
1  299.80000000
2  649.74000000
3  922.61000000
4  808.32540000
5  562.81362000
6  410.07991000
7  358.91942420
8  331.78020406
9  288.53786437
10 237.05769994
11 193.99105150
12 163.12702831
13 139.42906867
14 118.54110378
15  99.66403258
16  83.47154617
17  70.12395704
18  59.09140976
19  49.76856833
20  41.80868272
21  35.04996990
22  29.35710615
23  24.56986509
24  20.53155703
25  17.11645582
26  14.22899989
27  11.79098966
28   9.73376339
29   7.99726002
30   6.53075002
31   5.29213163
32   4.24619411
33   3.36309994
34   2.61748577
35   1.98789578
36   1.45625502
37   1.00733439
38   0.62827472
39   0.30820619
40   0.03794508
Show code
# Exibindo os resultados
  resultados <- data.frame(Periodo = 1:40, x, r, m, nv, y)

# Criando uma string
  y_values_list <- paste(resultados$y, collapse = ", ")

# Exibindo a string
  print(y_values_list)
[1] "299.8, 649.74, 922.61, 808.3254, 562.81362, 410.07991, 358.9194242, 331.78020406, 288.53786437, 237.0576999366, 193.99105150178, 163.12702831023, 139.429068672962, 118.541103775814, 99.6640325836356, 83.4715461665124, 70.12395703732, 59.0914097603875, 49.7685683308407, 41.8086827186021, 35.0499698974769, 29.3571061542066, 24.5698650919163, 20.5315570346214, 17.1164558232517, 14.2289998921823, 11.7909896610171, 9.73376338608289, 7.99726002040977, 6.5307500191365, 5.29213162559909, 4.24619410524074, 3.36309993961058, 2.61748577129534, 1.98789577894433, 1.45625502375225, 1.00733438584773, 0.628274717708676, 0.308206194619637, 0.0379450776554265"
Show code
# Convertendo em numérica
  y_values_numeric <- as.numeric(unlist(strsplit(y_values_list, ", ")))

# Condições iniciais
  initial_conditions <- c(y = y[1])
  
# Fazendo epidem.flowField
  epidem.flowField <- flowField(epidem,
                                xlim = c(0, 40),
                                ylim = c(0, 1000),
                                parameters = list(tr = tr, nv = nv0, x = x0),
                                points = 30,
                                system = "one.dim", add = FALSE, xlab = "t")

# Simulação de trajetória

    # Iterações
      times <- seq(0, 40, by = 1)
      
    # Trajetória levando em conta as iterações
      epidem.trajectory <- function(times, parameters, initial_conditions) {
      ode_result <- ode(y = initial_conditions,
                        times = times,
                        func = epidem,
                        parms = parameters) 
  
      return(as.data.frame(ode_result))
}

# Adicionando a trajetória ao gráfico usando a função lines
  for (i in 2:40) {
    ode_result <- epidem.trajectory(times = c(i-1, i),
                                    parameters = list(tr = tr, nv = nv[i-1], 
                                    x = x[i-1]), 
                                    initial_conditions = c(y = y[i-1]))
    lines(ode_result$time, ode_result$y, col = "red", type = "h")
}
Show code
  epidem.stability <- stability(epidem,
                                ystar = 0, 
                                parameters = list(tr = tr, nv = nv0, x = x0), 
                                system = "one.dim")
    #Apesar da classificação, consegue-se dizer, por meios teóricos, que y tende a 0