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

## Vector Fields and Trajectories!

\n", "

In this notebook, we're going to learn how to make vector fields and combine them with trajectory curves.  We will use the rabbits and foxes example from Section 8.1 in Calculus in Context to illustrate this process. Notice the following things about the plot_vector_field command:

\n", "

The first argument is a list of derivatives (in this case, (Rprime, Fprime)).

\n", "

The next two arguments specify the x-range and y-range.

\n", "

We have to define the variables first, or else Sage doesn't know that R and F are variables.

\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": [ "R,F = var('R F')\n", "plot_vector_field( (0.1*R*(1 - R/10000) - 0.005*R*F, \n", " 0.00004*R*F - 0.04*F), (R,0,3000), (F,0,30) )" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Notice that you can do this with pretty much any derivatives (try it!).

\n", "

Now, we already know how to plot a trajectory---we just do it by using a modification of our Euler's Method code to rabbits and foxes, and instead of doing t vs. R and t vs. F, we do R vs. F.  Like so:

" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starttime = 0\n", "stoptime = 300\n", "nsteps = 3000\n", "R = 1000\n", "F = 7\n", "deltat = float((stoptime-starttime)/nsteps)\n", "R_sol = [R]\n", "F_sol = [F]\n", "tvals = [starttime]\n", "for step in range(nsteps):\n", " R_prime = 0.1*R_sol[-1]*(1 - R_sol[-1]/10000) - 0.005*R_sol[-1]*F_sol[-1]\n", " F_prime = 0.00004*R_sol[-1]*F_sol[-1] - 0.04*F_sol[-1]\n", " R_sol.append(R_sol[-1] + deltat*R_prime)\n", " F_sol.append(F_sol[-1] + deltat*F_prime)\n", " tvals.append(tvals[-1] + deltat)\n", "\n", "list_plot(zip(R_sol,F_sol), color='purple')" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

It's hard to figure out the values at various points, so we can put a grid on our image.

" ] }, { "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": [ "list_plot(zip(R_sol,F_sol), color='purple').show(gridlines=True)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Hm, it would be more useful if we could have a finer grid...

" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list_plot(zip(R_sol,F_sol), color='purple').show(gridlines=\"minor\")" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

(Also, instead of \"minor\" you can use things like [\"major\",None] to get lines in one direction but not the other.)

\n", "

Oh, wanna show those two plots on the same axes?  Let's go for it.

" ] }, { "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": [ "R,F = var('R F')\n", "show(plot_vector_field( (0.1*R*(1 - R/10000) - 0.005*R*F, \n", " 0.00004*R*F - 0.04*F), (R,0,3000), (F,0,30) ) + list_plot(zip(R_sol,F_sol), color='purple'))" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Notice that we had to re-declare R and F as variables because they'd just been defined as constants when we ran our Euler's Method code.

\n", "

Here's a slightly easier-to-read way of doing the same thing.

" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "R,F = var('R F')\n", "VFplot = plot_vector_field( (0.1*R*(1 - R/10000) - 0.005*R*F, \n", " 0.00004*R*F - 0.04*F), (R,0,3000), (F,0,30) ) \n", "RFplot = list_plot(zip(R_sol,F_sol), color='purple')\n", "AllPlot = VFplot + RFplot\n", "AllPlot.show(gridlines=\"minor\")" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

We can do even better, though:  We can show the evolution of the trajectory through the vector field over time.

" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@interact\n", "def rabbits_foxes_fixed(stoptime = (70,(10, 1200))):\n", " starttime = 0\n", " nsteps = 3000\n", " deltat = float((stoptime-starttime)/nsteps)\n", " rab_sol = [1000]\n", " fox_sol = [7]\n", " tvals = [starttime]\n", " for step in range(nsteps):\n", " rab_prime = .1*rab_sol[-1]*(1 - rab_sol[-1]/10000) - .005*rab_sol[-1]*fox_sol[-1]\n", " fox_prime = 0.00004*rab_sol[-1]*fox_sol[-1] - 0.04*fox_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", " rabfoxplot = list_plot(zip(rab_sol, fox_sol),plotjoined=True,color='blue', title='Rabbits vs. Foxes')\n", " R,F = var('R F')\n", " vecfieldplot = plot_vector_field( (0.1*R*(1 - R/10000) - 0.005*R*F, 0.00004*R*F - 0.04*F), (R,0,max(rab_sol)), (F,0,max(fox_sol)) )\n", " RFplot = rabfoxplot + vecfieldplot \n", " RFplot.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Better yet, we can also fiddle with the initial conditions:

" ] }, { "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": [ "@interact\n", "def rabbits_foxes_init(stoptime = (70,(10, 1200)),rab = input_box(1000, label = 'initial number of rabbits: '), fox = input_box(5, label = 'initial number of foxes: ')):\n", " starttime = 0\n", " nsteps = 3000\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 = .1*rab_sol[-1]*(1 - rab_sol[-1]/10000) - .005*rab_sol[-1]*fox_sol[-1]\n", " fox_prime = 0.00004*rab_sol[-1]*fox_sol[-1] - 0.04*fox_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", " rabfoxplot = list_plot(zip(rab_sol, fox_sol),plotjoined=True,color='blue', title='Rabbits vs. Foxes')\n", " R,F = var('R F')\n", " vecfieldplot = plot_vector_field( (0.1*R*(1 - R/10000) - 0.005*R*F, 0.00004*R*F - 0.04*F), (R,0,max(rab_sol)), (F,0,max(fox_sol)) )\n", " RFplot = rabfoxplot + vecfieldplot \n", " RFplot.show()" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

Still not enough for you?  Oh, ho, ho, let's fiddle with the parameters!

" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def rabbits_foxes(stoptime,nsteps,rab,fox,a,b,c,d,e):\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]\n", " fox_prime = d*rab_sol[-1]*fox_sol[-1] - e*fox_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", " rabfoxplot = list_plot(zip(rab_sol, fox_sol),plotjoined=True,color='blue', title='Rabbits vs. Foxes')\n", " R,F = var('R F')\n", " vecfieldplot = plot_vector_field( (a*R*(1 - R/b) - c*R*F, d*R*F - e*F), (R,0,max(rab_sol)), (F,0,max(fox_sol)) )\n", " RFplot = rabfoxplot + vecfieldplot \n", " RFplot.show()" ] }, { "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": [ "rabbits_foxes(300, 3000, 1000, 7, .1, 10000, .005, .00004, .04)" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "

You know what you can do next.... interact over time, over parameters... go for it.

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