04-08-2020, 08:13 AM
Chances are you're reading this from your "shelter in place" Covid-19 confinement, and thus will resonate with the concepts involved. This short program is a direct application of Jean-Marc Baillard's ODE routine "RK4C" from the Differential Equations module, available at TOS and the CL Library (id$ = "DIFF").
The SIR model calculates the values of Susceptible, Infected and Removed groups of population in a total population of N individuals. N = S + I + R.
The virus is modeled with two parameters that indicate the infection rate "a" and the recovery rate "b". These are the crux of the model, as they need to be expressed in the same units used by the ODEs to show the individual results - i.e. a=2.3 infected people per person and DAY, and b = 0.3 people recovered per DAY if we want to look at daily numbers.
Supposedly your local Covid-19 statistics could be used to estimate the parameters, but this is tricky since the reported infection cases are much lower than the actual ones. Besides, the Removed section includes both the recovered (cured) and the dead patients.
The ODEs are normalized by the total population size N, but this is transparent to the user and it's done by the routine itself. The ODE's are implemented in the LBL "d/dt" subroutine:
dS/dt = -a. I(t).S(t)
dI/dt = a. S(t). I(t) - b.I(t)
dR/dt = b.I(t)
So here you have it, play with your a,b parameters (XEQ A) and define your initial conditions N, So, Io (XEQ B) to start getting the results for the S,I,R sectors of population. The model assumes Ro=0. Remember that N=S+I+R at all times.
Program Listing:
Note: ARCLI is in the AMC_OS/X Module. It appends the integer part o X to ALPHA.
PS. Looks that the formatting codes don't work within a CODE block?
PSS. - SIR material abounds on the Web, here's one video that was useful to me to grasp some details:
https://www.youtube.com/watch?v=k6nLfCbAzgo&t=867s
The SIR model calculates the values of Susceptible, Infected and Removed groups of population in a total population of N individuals. N = S + I + R.
The virus is modeled with two parameters that indicate the infection rate "a" and the recovery rate "b". These are the crux of the model, as they need to be expressed in the same units used by the ODEs to show the individual results - i.e. a=2.3 infected people per person and DAY, and b = 0.3 people recovered per DAY if we want to look at daily numbers.
Supposedly your local Covid-19 statistics could be used to estimate the parameters, but this is tricky since the reported infection cases are much lower than the actual ones. Besides, the Removed section includes both the recovered (cured) and the dead patients.
The ODEs are normalized by the total population size N, but this is transparent to the user and it's done by the routine itself. The ODE's are implemented in the LBL "d/dt" subroutine:
dS/dt = -a. I(t).S(t)
dI/dt = a. S(t). I(t) - b.I(t)
dR/dt = b.I(t)
So here you have it, play with your a,b parameters (XEQ A) and define your initial conditions N, So, Io (XEQ B) to start getting the results for the S,I,R sectors of population. The model assumes Ro=0. Remember that N=S+I+R at all times.
Program Listing:
Code:
9:41AM 04/08
01*[b]LBL "SIR"[/b]
02*[u]LBL A[/u]
03 RCL 12
04 "a="
05 ARCL X
06 >"?"
07 PROMPT
08 STO 12
09 RCL 13
10 "b="
11 ARCL X
12 >"?"
13 PROMPT
14 STO 13
15*[u]LBL B[/u]
16 RCL 14
17 "N="
18 ARCL X
19 >"?"
20 PROMPT
21 STO 14
22 RCL 02
23 *
24 "S0="
25 ARCL X
26 >"?"
27 PROMPT
28 RCL 14
29 /
30 STO 02
31 RCL 03
32 RCL 14
33 *
34 "I0="
35 ARCL X
36 >"?"
37 PROMPT
38 RCL 14
39 /
40 STO 03
41 0
42 STO 04
43 STO 01
44 ,1
45 STO 05
46 10
47 STO 06
48*LBL C
49 RCL 04
50 RCL 14
51 *
52 RCL 03
53 RCL 14
54 *
55 RCL 02
56 RCL 14
57 *
58 RCL 01
59 "S"
60 ARCLI
61 "`="
62 ARCL Y
63 PROMPT
64 "I"
65 ARCLI
66 >"="
67 ARCL Z
68 PROMPT
69 "R"
70 ARCLI
71 >"="
72 ARCL T
73 PROMPT
74 "d/dT"
75 ASTO 00
76 [color=#0000CD][b]XROM "RK4C"[/b][/color]
77 GTO C
78 RTN
79*[b]LBL "d/dT"[/b]
80 RDN
81 X<>Y
82 *
83 LASTX
84 RCL 13
85 *
86 X<>Y
87 RCL 12
88 *
89 ENTER^
90 CHS
91 X<>Y
92 RCL Z
93 ST- Y
94 END
Note: ARCLI is in the AMC_OS/X Module. It appends the integer part o X to ALPHA.
PS. Looks that the formatting codes don't work within a CODE block?
PSS. - SIR material abounds on the Web, here's one video that was useful to me to grasp some details:
https://www.youtube.com/watch?v=k6nLfCbAzgo&t=867s