{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

## There and Back

\n", "

Suppose you want to measure the accumulated error from doing Euler's method (remember, that's the official name for our S-I-R calculations).  What you want to have the computer do is

\n", "

(1) run the S-I-R code from starttime to stoptime

\n", "

(2) note the ending S-I-R from that procedure

\n", "

(3) use them as the initial values for running the S-i-R code from old stoptime to old starttime

\n", "

Here is code to do step (1), then code to do step (3):

" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starttime = 0\n", "stoptime = 1\n", "nsteps = 10\n", "S = 45400\n", "I = 2100\n", "R = 2500\n", "t = starttime\n", "deltat = float((stoptime-starttime)/nsteps)\n", "for step in range(nsteps):\n", " S_prime = -.00001*S*I\n", " I_prime = .00001*S*I - I/14\n", " R_prime = I/14\n", " deltaS = S_prime*deltat\n", " deltaI = I_prime*deltat\n", " deltaR = R_prime*deltat\n", " t = t + deltat\n", " S = S + deltaS\n", " I = I + deltaI\n", " R = R + deltaR\n", "print(t,S,I,R)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starttime = 1\n", "stoptime = 0\n", "nsteps = 10\n", "S = (put the output value for S here)\n", "I = (put the output value for I here)\n", "R = (put the output value for R here)\n", "t = starttime\n", "deltat = float((stoptime-starttime)/nsteps)\n", "for step in range(nsteps):\n", " S_prime = -.00001*S*I\n", " I_prime = .00001*S*I - I/14\n", " R_prime = I/14\n", " deltaS = S_prime*deltat\n", " deltaI = I_prime*deltat\n", " deltaR = R_prime*deltat\n", " t = t + deltat\n", " S = S + deltaS\n", " I = I + deltaI\n", " R = R + deltaR\n", "print(t,S,I,R)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Oh, but there's a more efficient way to do it.  The thing is that you can run it without re-setting S, I, and R!  See below.

" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starttime = 0\n", "stoptime = 1\n", "nsteps = 10\n", "S = 45400\n", "I = 2100\n", "R = 2500\n", "t = starttime\n", "deltat = float((stoptime-starttime)/nsteps)\n", "for step in range(nsteps):\n", " S_prime = -.00001*S*I\n", " I_prime = .00001*S*I - I/14\n", " R_prime = I/14\n", " deltaS = S_prime*deltat\n", " deltaI = I_prime*deltat\n", " deltaR = R_prime*deltat\n", " t = t + deltat\n", " S = S + deltaS\n", " I = I + deltaI\n", " R = R + deltaR\n", "starttime = 1\n", "stoptime = 0\n", "nsteps = 10\n", "t = starttime\n", "deltat = float((stoptime-starttime)/nsteps)\n", "for step in range(nsteps):\n", " S_prime = -.00001*S*I\n", " I_prime = .00001*S*I - I/14\n", " R_prime = I/14\n", " deltaS = S_prime*deltat\n", " deltaI = I_prime*deltat\n", " deltaR = R_prime*deltat\n", " t = t + deltat\n", " S = S + deltaS\n", " I = I + deltaI\n", " R = R + deltaR\n", "print(t,S,I,R)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

## Basic Data Plotting

\n", "

In order to plot data, we have to actually have data to plot---which, for a computer, means that all the data has to be in one place. Thus, we make lists of values as we go. In the code, this happens by

\n", "
\n", "
• Creating lists for the time data, S data, I, data, and R data (written in square brackets, i.e. [ ])
• \n", "
• For each time value, use the last value in the list (the [-1] value) to compute a new value, and append the new value to the list.
• \n", "
\n", "
This code has no output---it just loads data into lists.
\n", "

" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starttime = 0\n", "stoptime = 3\n", "nsteps = 6\n", "S = 45400\n", "I = 2100\n", "R = 2500\n", "deltat = float((stoptime-starttime)/nsteps)\n", "S_sol = [S]\n", "I_sol = [I]\n", "R_sol = [R]\n", "tvals = [starttime]\n", "for step in range(nsteps):\n", " S_prime = -.00001*S_sol[-1]*I_sol[-1]\n", " I_prime = .00001*S_sol[-1]*I_sol[-1] - I_sol[-1]/14\n", " R_prime = I_sol[-1]/14\n", " S_sol.append(S_sol[-1] + deltat*S_prime)\n", " I_sol.append(I_sol[-1] + deltat*I_prime)\n", " R_sol.append(R_sol[-1] + deltat*R_prime)\n", " tvals.append(tvals[-1] + deltat)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Now we want to plot the data.  First we will plot the S data, with time on the horizontal axis and number of people on the vertical axis.  The zip command makes coordinate pairs, with the x-value from the first list and the y-value from the second list.

" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(zip(tvals,S_sol))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

If we want to, we can connect the dots.

" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(zip(tvals,S_sol),plotjoined=True)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

We can also label the graph.

" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(zip(tvals,S_sol),plotjoined=True,title='Time in days vs. Number of Susceptibles')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

You can alter the code above to see other plots---time vs. infecteds, susceptibles vs. recovereds, etc.

\n", "

For fancier plotting (multiple graphs at the same time, for example), please see the next notebook!

" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 8.0", "language": "", "name": "sagemath" }, "language": "python", "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }