Sir Backwards

R(t) = 1 / 5days I(t)

R2 = R1 + h 1/5daysI2
S2 = S1 + h(- 5 * 10^-9/ day*Person I1,S1)

Implicit solvers: use it or not?

Trapezoidal Stability
Cases – Moderate value of h, hk/2 = 0.1
hk/2 = 1
High value of h, e.g. hk/2 = g

Models with flows between different compartments

from xxx import *

h = 0.5
end_time = 60.
num_steps = int(end_time / h)
times = h * numpy.array(range(num_steps + 1))

def waning():
	transmission_coeff = 5e-9
	infectious_time = 5.

	waning_time

	s = numpy.zeros(num_steps + 1)
	i = numpy.zeros(num_steps + 1)
	r = numpy.zeros(num_steps + 1)

	s[0] = 1e8 - 1e6 - 1e5
	i[0] = 1e5
	r[0] = 1e6

	for step in range(num_steps):
		s2i = h * transmission_coeff * s[step] * i[step]
		i2r = h / infectious_time * i[step]

		s[step + 1] = s[step] - s2i
		i[step + 1] = i[step] + s2i - i2r
		r[step + 1] = r[step] + i2r

		return s, i, r

s, i, r = waning()

@show_plot
def plot_me():
    s_plot = matplotlib.pyplot.plot(times, s, label = 'S')
    i_plot = matplotlib.pyplot.plot(times, i, label = 'I')
    r_plot = matplotlib.pyplot.plot(times, r, label = 'R')
    matplotlib.pyplot.legend(('S', 'I', 'R'), loc = 'upper right')

    axes = matplotlib.pyplot.gca()
    axes.set_xlabel('Time in days')
    axes.set_ylabel('Number of persons')
    matplotlib.pyplot.xlim(xmin = 0.)
    matplotlib.pyplot.ylim(ymin = 0.)
plot_me()