Question about Euler's Method and the SIR epidemic model using a spreadsheet

I am creating a spreadsheet to figure out how to correctly use Euler's Method with the SIR Model for Spread of Disease. I took a random question from the text Calculus in Context by Callahan and Hoffman.

Given that $S(t)$ is the number of susceptible individuals, $I(t)$ is the number of infected individuals, and $R(t)$ is the number of recovered individuals, I am using the following differential equations:

$\frac{dS}{dt} = -bSI$

$\frac{dI}{dt} = bSI-kI$

$\frac{dR}{dt} = kI$

So, with Euler's method, I used the following idea to get each successive iteration:

$S_{new}=S-(bSI)\Delta t$

$I_{new}=I+(bSI-kI)\Delta t$

$R_{new}=R+(kI)\Delta t$

The particular problem I am working with has the following initial condition:

$S(0)=35,400$

$I(0)=13,500$

$R(0)=22,100$

Also, the text set $b = 0.00001$ and $k=0.08$. I made $\Delta t = 1$ to calculate $S$, $I$, and $R$ each day.

Here is a screenshot of the spreadsheet showing my results for the first 10 days:

Screenshot of the first 10 days

Here is what I have entered in the cells

Note that $b$ is stored as 0.00001, $k$ is stored as 0.08, and $t$ is the $\Delta t$, which is stored as 1.

In cell B3 I have entered: =B2-b*B2*C2*t

In cell C3, I have entered: =C2+(b*B2*C2-k*C2)*t

In cell D3, I have entered: D2+k*C2*t

For the rest of the cells up to 40 days, I just dragged down the same formulas. So, for B4, C4, and D4, for example, it would be the identical formulas I listed above except with B3, C3 and D3 inside.

• I am first making sure I have the right idea and am not totally missing something crucial.
• I am also very interested in relating this to the COVID-19. Correct me if I'm wrong, but I believe $k$ is 1/(average period of infectiousness). I am still confused about how to estimate a $b$ value. I'd like to relate this to the US population, but, if that is too big, maybe another specific state or city. If someone could help me a bit with this piece, that would be great!

Solution 1:

The Euler method with a larger step size $Δt$ is not really faithful to the original model $\dot u=f(u)$, it is closer to a solution of $$ \dot u=f_{Δt}(u)=f(u)-\frac{Δt}2f'(u)f(u) $$ plus higher order/degree terms. This means samples of the exact solution $u_{Δt}$ of this equation are $O(Δt^2)$ close to the points obtained from the Euler method, as the Taylor expansion of a time step shows, \begin{align} u_{Δt}(t+Δt)&=u_{Δt}(t)+f_{Δt}(u_{Δt}(t))Δt+\frac12f_{Δt}'f_{Δt}Δt^2+...\\ &=u_{Δt}(t)+[f(u_{Δt})-\frac12f'(u_{Δt})f(u_{Δt})Δt]Δt+\frac12f'(u_{Δt})f(u_{Δt})Δt^2+O(Δt^3) \\ &=u_{Δt}(t)+f(u_{Δt}(t))Δt+O(Δt^3). \end{align} This means that the time step is an Euler step plus an error $O(Δt^3)$ which accumulates over time to a global error $O(Δt^2)$.

dotted - original equation, lines - modified equation, crosses - Euler steps

The first order modification will contain cubic terms, higher order modification give higher degree terms, giving it a different quality relative to the original equation. The difference of the Euler-discretized numerical values and the original ODE can not simply be captured by an adaptation of the constants

So it would be more honest to formulate the simulation directly as discrete dynamical model, or to make the ODE integration more accurate, meaning either a reduced step size in the Euler method, 10 to 50 steps per day, or a higher order method with 1 to 4 steps per day. In Excel, this would require using VBA primarily, using the spreadsheet itself only to hold the input and output data.

To compare against, I produced the first 10 values at integer times for these different methods, first with time step 1 as in your computation, then Euler with a smaller step size and RK4 as the much more accurate method.

Euler with step size dt=1
-------------------------
[[35400.         13500.         22100.        ]
 [30621.         17199.         23180.        ]
 [25354.49421    21089.58579    24555.92      ]
 [20007.33640196 24749.57673484 26243.0868632 ]
 [15055.60532656 27721.34167145 28223.05300199]
 [10881.98953328 29677.25013102 30440.7603357 ]
 [ 7652.51428026 30532.54537356 32814.94034618]
 [ 5316.00688542 30426.44913851 35257.54397607]
 [ 3698.53475423 29609.80533862 37691.65990715]
 [ 2603.40581312 28336.14985264 40060.44433424]
 [ 1865.70084064 26806.96283691 42327.33632245]]

Euler with step size dt=0.1, every 10th value printed
-------------------------
[[35400.         13500.         22100.        ]
 [30372.89492612 17311.1994966  23315.90557728]
 [25061.76732461 21099.40208755 24838.83058784]
 [19930.44031117 24418.83970897 26650.71997986]
 [15373.94595865 26925.48421687 28700.56982448]
 [11608.55630897 28474.70354729 30916.74014374]
 [ 8664.12109087 29112.62797744 33223.25093169]
 [ 6447.61114318 28999.87376973 35552.51508708]
 [ 4816.88935953 28331.62601808 37851.48462239]
 [ 3630.2581671  27287.36046897 40082.38136393]
 [ 2768.74774096 26010.70062585 42220.55163319]]

RK4 with step size dt=0.25, every 4th value printed
-------------------------
[[35400.         13500.         22100.        ]
 [30348.86304515 17319.5147641  23331.62219075]
 [25040.13784695 21090.00792405 24869.85422901]
 [19937.63568967 24369.54314768 26692.82116265]
 [15422.7274815  26830.31557838 28746.95694012]
 [11695.20776528 28344.53042335 30960.26181137]
 [ 8775.26846677 28966.5348119  33258.19672133]
 [ 6568.54882662 28856.09841697 35575.35275641]
 [ 4936.33106926 28202.95105006 37860.71788068]
 [ 3741.40309106 27180.58529638 40078.01161256]
 [ 2868.45607877 25928.03780371 42203.50611751]]