Introduction
This section develops a Python program that solves the SI differential equation. Examining the output shows that the differential equation solutions agrees well with the actual solution. I call your attention to the passing of two parameters through the parameter list, as well as some formatting get the columns nice and neat. As we progress we will add more to this code.
In order to handle the restricted column width, use the slide bar at the bottom of the Python Code and output.
Python Code for the SI differential equation.
# ************** SI Model differential equation Solver
import numpy as np
from scipy.integrate import odeint
def f(y, t, parm): # returns the derivative of f(t) (y)
I = y[0] # get the current value of y
plam,Npop = parm # get the parameter
derivative = [plam*I*(1-I/Npop)] # calculate the derrivative
# print(I,plam,derivative) # if you want to see where odeint is looking
return derivative
plam = 0.2137 # growth rate at start (infections per day per infector)
Npop = 50000.0 # Total "population"
I0 = 1619.0 # Initial number of infected individuals
e=2.7182818284 # Euler's number (approx)
parm = [plam,Npop] #ODE solver wants a list of Parameters
y0 = [I0] #ODE also wants a list of initial values
# Make an array of the times at which we want a solution
Finaltime = 60.
Increment = 1.0
times = np.arange(0.0, Finaltime, Increment)
# Call the ODE solver
Solutions = odeint(f, y0, times, args=(parm,))
print(" day, Susceptible, infectious, formula, error ")
for i in range(len(times)): # printout and compare with actual solution
t = times[i]
I = Solutions[i,0]
S = Npop - I
val = Npop/(1+(Npop-I0)/I0 * e**(-1.0*plam*t) ) # the actual solution
error = (val-I)/val
print("{:4.1f}".format(t), "{:11.1f}".format(S),"{:11.1f}".format(I),
"{:11.1f}".format(val)," ","{:8.4E}".format(error))
Python Output
day, Susceptible, infectious, formula, error
0.0 48381.0 1619.0 1619.0 0.0000E+00
1.0 48010.6 1989.4 1989.4 9.7049E-09
2.0 47559.8 2440.2 2440.2 1.8578E-08
3.0 47013.1 2986.9 2986.9 3.4326E-08
4.0 46353.4 3646.6 3646.6 5.0050E-08
5.0 45561.7 4438.3 4438.3 6.6686E-08
6.0 44618.1 5381.9 5381.9 5.8821E-08
7.0 43502.5 6497.5 6497.5 5.7928E-08
8.0 42196.1 7803.9 7803.9 5.4188E-08
9.0 40683.2 9316.8 9316.8 2.8894E-08
10.0 38953.9 11046.1 11046.1 1.9692E-08
11.0 37006.1 12993.9 12993.9 5.7516E-09
12.0 34848.4 15151.6 15151.6 -6.9563E-09
13.0 32501.8 17498.2 17498.2 -1.5389E-08
14.0 30000.4 19999.6 19999.6 -8.2386E-09
15.0 27390.2 22609.8 22609.8 2.0009E-09
16.0 24726.2 25273.8 25273.8 6.6028E-09
17.0 22068.5 27931.5 27931.5 8.5012E-09
18.0 19476.3 30523.7 30523.7 5.4730E-09
19.0 17003.3 32996.7 32996.7 -3.3017E-09
20.0 14693.1 35306.9 35306.9 -1.4317E-08
21.0 12577.1 37422.9 37422.9 -2.3828E-08
22.0 10673.8 39326.2 39326.2 -2.4812E-08
23.0 8989.3 41010.7 41010.7 -2.2251E-08
24.0 7519.8 42480.2 42480.2 -1.6971E-08
25.0 6253.9 43746.1 43746.1 -1.0817E-08
26.0 5175.1 44824.9 44824.9 -6.1004E-09
27.0 4264.3 45735.7 45735.7 -3.1126E-09
28.0 3501.3 46498.7 46498.7 -1.1449E-09
29.0 2866.2 47133.8 47133.8 -5.8592E-10
30.0 2340.5 47659.5 47659.5 6.0688E-11
31.0 1907.4 48092.6 48092.6 -4.0188E-10
32.0 1551.8 48448.2 48448.2 -1.6644E-09
33.0 1260.7 48739.3 48739.3 -2.4823E-09
34.0 1023.1 48976.9 48976.9 -3.0223E-09
35.0 829.5 49170.5 49170.5 -3.0706E-09
36.0 672.1 49327.9 49327.9 -3.0638E-09
37.0 544.2 49455.8 49455.8 -3.0746E-09
38.0 440.4 49559.6 49559.6 -3.1083E-09
39.0 356.2 49643.8 49643.8 -2.6492E-09
40.0 288.1 49711.9 49711.9 -2.3356E-09
41.0 232.9 49767.1 49767.1 -1.8957E-09
42.0 188.3 49811.7 49811.7 -1.5676E-09
43.0 152.2 49847.8 49847.8 -1.2665E-09
44.0 123.0 49877.0 49877.0 -9.8527E-10
45.0 99.3 49900.7 49900.7 -7.8722E-10
46.0 80.3 49919.7 49919.7 -7.7493E-10
47.0 64.8 49935.2 49935.2 -1.0414E-09
48.0 52.4 49947.6 49947.6 -1.5342E-09
49.0 42.3 49957.7 49957.7 -1.7176E-09
50.0 34.2 49965.8 49965.8 -1.7072E-09
51.0 27.6 49972.4 49972.4 -1.5655E-09
52.0 22.3 49977.7 49977.7 -1.3366E-09
53.0 18.0 49982.0 49982.0 -1.1382E-09
54.0 14.5 49985.5 49985.5 -1.0085E-09
55.0 11.7 49988.3 49988.3 -9.1308E-10
56.0 9.5 49990.5 49990.5 -8.1465E-10
57.0 7.7 49992.3 49992.3 -7.2737E-10
58.0 6.2 49993.8 49993.8 -6.1825E-10
59.0 5.0 49995.0 49995.0 -5.1561E-10
Next steps
There is another way of implementing the SI model on a computer. Here is a link takes to a Monte Carlo technique. https://davidalarrabee.com/?page_id=873