Euler Estimation and SIR Model

Arrays can be used to implement more than just matrices. We can also implement other mathematical procedures like Euler estimation.

Single-Variable Euler Estimator

To start, build a single-variable Euler estimation class as follows:

>>> def derivative(t):
        return t+1
>>> euler = EulerEstimator(derivative)

>>> initial_point = (1,4)
>>> euler.eval_derivative(initial_point) # evaluates derivative at point (1,4)

>>> step_size = 0.5
>>> num_steps = 4
>>> euler.estimate_points(initial_point, step_size, num_steps)
    (1, 4),      # starting point
    (1.5, 5),    # after 1st step
    (2, 6.25),   # after 2nd step
    (2.5, 7.75), # after 3rd step
    (3, 9.5)     # after 4th step

Then, use your Euler estimator to plot several solution curves to the following differential equation on the interval $x \in [0, 5].$ (Your Euler estimator generates a list of points, and then you can use that list of points to generate a plot.)

$\begin{align*} \dfrac{\textrm dy}{\textrm dx} = x - 2 \end{align*}$

For one curve, use the initial condition $y(0)=-2.$ For another curve, use $y(0)=-1.$ Then another curve with $y(0)=0,$ another with $y(0) = 1,$ and another with $y(0)=2.$ All $5$ of these curves can go on the same plot, and based on your knowledge of calculus, you should be able to tell if your plots look right.

Multivariable Euler Estimator

Once you’ve impelmented a single-variable Euler estimator, you can generalize it to simulate systems of differential equations. For example, consider the following system:

$\begin{align*} a'(t) &= a(t) + 1 \\ b'(t) &= a(t) + b(t) \\ c'(t) &= 2b(t) + 3t \end{align*}$

To simulate this system starting with the initial state $a(-0.4) = -0.45,$ $b(-0.4) = -0.05,$ $c(-0.4) = 0,$ construct a multivariable Euler estimator as follows:

>>> initial_state = {'a': -0.45, 'b': -0.05, 'c': 0}

>>> initial_point = (-0.4, initial_state) # points take form (t, state)

>>> def da_dt(t, state):
        return state['a'] + 1

>>> def db_dt(t, state):
        return state['a'] + state['b']
>>> def dc_dt(t, state):
        return 2 * state['b'] + 3 * t
>>> derivatives = {
    'a': da_dt,
    'b': db_dt,
    'c': dc_dt
>>> euler = EulerEstimator(derivatives)

>>> euler.eval_derivative_at_point(initial_point)
{'a': 0.55, 'b': -0.5, 'c': -1.3}

>>> step_size = 2
>>> num_steps = 3
>>> euler.estimate_points(initial_point, step_size, num_steps)
   (-0.4, {'a': -0.45, 'b': -0.05, 'c': 0}),  # starting point 
   (1.6, {'a': 0.65, 'b': -1.05, 'c': -2.6)), # after 1st step
   (3.6, {'a': 3.95, 'b': -1.85, 'c': 2.8)),  # after 2nd step
   (5.6, {'a': 13.85, 'b': 2.35, 'c': 17))    # after 3rd step

SIR Model

One of the simplest ways to model the spread of disease using differential equations is the SIR model. Let’s construct a SIR model and use our multivariable Euler estimator to plot the solution curves.

The SIR model assumes three sub-populations: susceptible, infected, and recovered.

  • The number of susceptible people $(S)$ decreases at a rate proportional to the rate of meeting between susceptible and infected people (because susceptible people have a chance of catching the disease when they come in contact with infected people).
  • The number of infected people $(I)$ increases at a rate proportional to the rate of meeting between susceptible and infected people (because susceptible people become infected after catching the disease), and decreases at a rate proportional to the number of infected people (as the diseased people recover).
  • The number of recovered people $(R)$ increases at a rate proportional to the number of infected people (as the diseased people recover).

First, write a system of differential equations using the following assumptions:

  • There are initially $1000$ susceptible people and $1$ infected person.
  • The number of meetings between susceptible and infected people each day is proportional to the product of the numbers of susceptible and infected people, by a factor of $0.01.$ The transmission rate of the disease is $3\%.$ (In other words, $3\%$ of meetings result in transmission.)
  • Each day, $2\%$ of infected people recover.
$\begin{align*} \begin{cases} \dfrac{\textrm{d}S}{\textrm{d}t} &= \_\_\_, \quad S(0) = \_\_\_ \\ \dfrac{\textrm{d}I}{\textrm{d}t} &= \_\_\_, \quad I(0) = \_\_\_ \\ \dfrac{\textrm{d}R}{\textrm{d}t} &= \_\_\_, \quad R(0) = \_\_\_ \end{cases} \end{align*}$

If you’ve written the system correctly, then at $t=0$ you should have

$\begin{align*} \dfrac{\textrm{d}S}{\textrm{d}t} = -0.3, \quad \dfrac{\textrm{d}I}{\textrm{d}t} = 0.28, \quad \dfrac{\textrm{d}R}{\textrm{d}t} = 0.02 \, . \end{align*}$

Then use your multivariable Euler estimator to simulate the solution curve, and plot the result.

  • Be sure to choose the step size small enough that the simulation doesn't blow up, but large enough that it doesn't take long to run.
  • Likewise, choose an interval that displays all the main features of the differential equation, i.e. an interval that shows the behavior of the curves until they start to asymtpote off.

If your plot is correct, then you should be able to easily describe what is happening in the plot and why it is happening.

Leave a Comment