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

## Let's Play with the May Model!

\n", "

Here is a function defined in Sage to run an Euler approximation to a solution of the May model for rabbit/fox populations, and to plot the resulting data.  You'll notice that it's a multivariable functionit takes 10 variable inputs (ending time, number of steps, initial numbers of rabbits and foxes, and the 6 parameters of the model)---and this makes it very flexible.  Evaluate the cell below in order to be able to use the MayModel function.

\n", "

" ] }, { "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": [ "def may_model_graphs(stoptime,nsteps,rab,fox,a,b,c,d,e,f):\n", " starttime = 0\n", " deltat = float((stoptime-starttime)/nsteps)\n", " rab_sol = [rab]\n", " fox_sol = [fox]\n", " tvals = [starttime]\n", " for step in range(nsteps):\n", " rab_prime = a*rab_sol[-1]*(1 - rab_sol[-1]/b) - c*rab_sol[-1]*fox_sol[-1]/(rab_sol[-1] + d)\n", " fox_prime = e*fox_sol[-1]*(1 - fox_sol[-1]/(f*rab_sol[-1]))\n", " rab_sol.append(rab_sol[-1] + deltat*rab_prime)\n", " fox_sol.append(fox_sol[-1] + deltat*fox_prime)\n", " tvals.append(tvals[-1] + deltat)\n", " rabplot = list_plot(zip(tvals,rab_sol),plotjoined=True,color='blue', title='Time in days vs. Rabbits(blue) and Foxes(red)')\n", " foxplot = list_plot(zip(tvals,fox_sol),plotjoined=True, color='red')\n", " Mayplot = rabplot + foxplot \n", " Mayplot.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Go ahead and try it out for a sample set of parameters, initial values, etc.  (The sample below corresponds to problem 6c.)

" ] }, { "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": [ "may_model_graphs(30, 1000, 20, 10, .6, 10, .5, 1, .1, 2)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Now, here's an interact that allows you to play with the value of the parameter c, a la problem 6e.  (You'll see why this is called an 'interact' when you run the code...)

" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@interact\n", "def may_model_c_int(c= input_box(0.5, label = 'parameter c: ')):\n", " starttime = 0\n", " stoptime = 200\n", " nsteps = 1500\n", " deltat = float((stoptime-starttime)/nsteps)\n", " rab_sol = [20]\n", " fox_sol = [20]\n", " tvals = [starttime]\n", " for step in range(nsteps):\n", " rab_prime = .6*rab_sol[-1]*(1 - rab_sol[-1]/10) - c*rab_sol[-1]*fox_sol[-1]/(rab_sol[-1] + 1)\n", " fox_prime = .1*fox_sol[-1]*(1 - fox_sol[-1]/(2*rab_sol[-1]))\n", " rab_sol.append(rab_sol[-1] + deltat*rab_prime)\n", " fox_sol.append(fox_sol[-1] + deltat*fox_prime)\n", " tvals.append(tvals[-1] + deltat)\n", " rabplot = list_plot(zip(tvals,rab_sol),plotjoined=True,color='blue', title='Time in days vs. Rabbits(blue) and Foxes(red)')\n", " foxplot = list_plot(zip(tvals,fox_sol),plotjoined=True, color='red')\n", " Mayplot = rabplot + foxplot \n", " Mayplot.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Why stop there?  If you want to alter any or all of the parameters, use the following interact...

" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@interact\n", "def may_model_graph_int(stoptime = input_box(30, label = 'ending time: '),nsteps = input_box(1000, label = 'number of steps: '),rab = input_box(20, label = 'initial number of centirabbits: '), fox = input_box(10, label = 'initial number of foxes: '),a = input_box(0.6, label = 'parameter a: '),b = input_box(10, label = 'parameter b: '),c = input_box(0.5, label = 'parameter c: '),d = input_box(1, label = 'parameter d: '),e = input_box(0.1, label = 'parameter e: '),f = input_box(2, label = 'parameter f: ')):\n", " starttime = 0\n", " deltat = float((stoptime-starttime)/nsteps)\n", " rab_sol = [rab]\n", " fox_sol = [fox]\n", " tvals = [starttime]\n", " for step in range(nsteps):\n", " rab_prime = a*rab_sol[-1]*(1 - rab_sol[-1]/b) - c*rab_sol[-1]*fox_sol[-1]/(rab_sol[-1] + d)\n", " fox_prime = e*fox_sol[-1]*(1 - fox_sol[-1]/(f*rab_sol[-1]))\n", " rab_sol.append(rab_sol[-1] + deltat*rab_prime)\n", " fox_sol.append(fox_sol[-1] + deltat*fox_prime)\n", " tvals.append(tvals[-1] + deltat)\n", " rabplot = list_plot(zip(tvals,rab_sol),plotjoined=True,color='blue', title='Time in days vs. Rabbits(blue) and Foxes(red)')\n", " foxplot = list_plot(zip(tvals,fox_sol),plotjoined=True, color='red')\n", " Mayplot = rabplot + foxplot \n", " Mayplot.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

And why even stop there?  You could modify any of the above code (defining the function, an interact that lets you vary one parameter, an interact that lets you vary lots of parameters) to work with the Lotka-Volterra model, or even the S-I-R model... Let your imagination run wild!

" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 10, "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 }