SIRSERIES

In review

This page runs a BASIC program that finds the solutions of the SIR model as power series,

\[ S(t) = s_0 + s_1 t + s_2 t^2 + \cdots, \qquad I(t) = i_0 + i_1 t + i_2 t^2 + \cdots, \]

instead of stepping forward with Euler’s method. The program is rewritten in Python and runs right here in your browser. Edit the code and press Run to run it again; Reset brings back the original program.

The model is the same one as in SIRPLOT: \(S' = -aSI\) and \(I' = aSI - bI\), with \(a = 0.00001\), \(b = 1/14\), \(S(0) = 45{,}400\), and \(I(0) = 2100\). Putting the series into these equations and matching the coefficients of each power of t gives a rule for each coefficient in terms of the ones before it. The coefficient of \(t^{k-1}\) in the product \(SI\) is \(s_0 i_{k-1} + s_1 i_{k-2} + \cdots + s_{k-1} i_0\), which is what the inner loop adds up in sum. Then

\[ s_k = -\frac{a \cdot \text{sum}}{k}, \qquad i_k = -s_k - \frac{b\, i_{k-1}}{k}. \]

The arrays s and i hold the coefficients; DIM s(50), i(50) makes room for 51 of each, numbered 0 to 50. In Python these are lists, s = [0] * 51, written with square brackets: s[k] instead of s(k).

The original program only computes the coefficients and doesn’t print anything, so the Python adds one line, marked in a comment, that prints k, \(s_k\), and \(i_k\). Notice how quickly they shrink: \(s_1 = -953.4\), \(s_2 \approx -172.4\), \(s_3 \approx -18.0\), and so on.

The extension. Check Include the extension to add the textbook’s extension to the end of the program. It defines SUS(x) and INF(x), which add up the series through the \(x^{50}\) term, and graphs both from 0 to 10: S in blue and I in red. Unchecking the box removes it again; your changes to the main program stay put. The extension is only an outline in the textbook (“Plot the line from … to …”, with the graphics setup left open), so the Python fills in line(...) commands and a window(...) that fits the numbers. The curves should look like the ones from SIRPLOT. Try widening the window and the loop to 20 or 30 to see what happens to a power series far from its center.

The textbook prints the first line as a = .00001q; the stray q is a typo, so the Python uses 0.00001.

Notes from the textbook

Add notes here.

 1  DEFDBL A-Z
 2  DIM s(50), i(50)
 3  a = .00001q
 4  b = 1 / 14
 5  s(0) = 45400
 6  i(0) = 2100
 7  FOR k = 1 TO 50
 8          sum = 0
 9          FOR j = 0 TO k - 1
10                  sum = sum + s(j) * i(k - j - 1)
11          NEXT j
12          s(k) = -a * sum / k
13          i(k) = -s(k) - b * i(k - 1) / k
14  NEXT k

The textbook gives this part in outline: “Plot the line from … to …” stands for a LINE command, and the graphics setup is left to you.

 1  DEF SUS(x)
 2       Sum = S(0)
 3       FOR j = 1 TO 50
 4           Sum = Sum + S(j) * x ^ j
 5       NEXT j
 6       SUS = Sum
 7  END DEF
 8  DEF INF(x)
 9       Sum = I(0)
10       FOR j = 1 TO 50
11           Sum = Sum + I(j) * x ^ j
12       NEXT j
13       INF = Sum
14  END DEF
15  FOR x = 0 TO 10 STEP .01
16       Plot the line from (x, SUS(x)) to (x + .01, SUS(x + .01))
17       Plot the line from (x, INF(x)) to (x + .01, INF(x + .01))
18  NEXT x

Output

Everything the program prints appears here, one line for each print(...), like BASIC's PRINT. Numbers printed on the same line are lined up in columns.

Loading Python. The first run takes a few seconds.

Python code

Press Run (or Ctrl/⌘ + Enter) to run the program again. Press Esc, then Tab, to leave the editor.

Loading Python…