Vector Direction Fields

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Special Mathematica Functions used in this NoteBook

Off[General::spell1];
BeginPackage["Graphics`ArrowedDirectionField`"];


shadeDirField::usage =
	"shadeDirField[f[u,v],g[u,v],{u,umin,umax,unum},
	{v,vmin,vmax,vnum}, (options) ]  
	plots a 2 dimn'l direction field determined by dx/dt=f 
	and dy/dt=g. unum and vnum determine number of arrows 
	plotted along x-axis and y-axis, resp. The shade of 
	the arrows show their relative length/strength.
	Options are passed to Show[]";

unitDirField::usage=
	"unitDirField[f[u,v],g[u,v],{u,umin,umax,unum},
	{v,vmin,vmax,vnum}, (options) ]  
	plots a 2 dimn'l direction field determined by dx/dt=f 
	and dy/dt=g. unum and vnum determine number of arrows 
	plotted along u-axis and v-axis, resp. 
	Options are passed to Show[]";
	
Begin["`Private`"];
	
Arrow[r_,theta_,pt_]:=
Module[{P1,P2,P},
	P1={6r/7 Cos[Pi/10 + theta] + pt[[1]], 
		6r/7 Sin[Pi/10 + theta] + pt[[2]]};
	P2={6r/7 Cos[-Pi/10 + theta]+ pt[[1]], 
		6r/7 Sin[-Pi/10 + theta]+ pt[[2]]};
	P={r Cos[theta]+ pt[[1]], r Sin[theta]+ pt[[2]]};
	{Polygon[{P1,P2,P}],Line[{pt,P}]} 
   ];
   
Theta[{x0_,y0_},{x1_,y1_}]:=
	If[x1 == x0 ,If[ y0 < y1, Pi/2, If[y0 > y1, 3Pi/2,0]],
	If[x1 < x0, Pi + ArcTan[(y1-y0)/(x1-x0)],
		ArcTan[(y1-y0)/(x1-x0)] ] ];
		
shadeDirField[f_,g_,{x_,xmin_,xmax_,xnum_},
	{y_,ymin_,ymax_,ynum_}, opts___]:=
Module[{points,values,linelength,M,m,i,j,dx,dy,xx},
		
gray[xx_]:=If[(M - xx)/(M-m) < .8, ((M - xx)/(M-m))^6,
	    (2.5((M - xx)/(M-m))-1.738)];
	   
	dx = N[(xmax - xmin)/(xnum-1)] ;
	dy = N[(ymax - ymin)/(ynum-1)] ;
	linelength = Min[dx, dy]/2 ;
	points=Flatten[Table[{i,j}, 
		{i,xmin,xmax,dx},{j,ymin,ymax,dy}], 1];
	values=Flatten[Table[{f[i,j],g[i,j]}, 
		{i,xmin,xmax,dx},{j,ymin,ymax,dy}], 1];
	norms=Table[N[Sqrt[values[[i]][[1]]^2 + 
		values[[i]][[2]]^2]],{i,1, Length[values]}] ;
	M = Max[norms] ; m = Min[norms] ;
	Show[Graphics[Table[{GrayLevel
		[gray[norms[[i]]]],
		Arrow[linelength, Theta[{0,0},values[[i]]], 
		points[[i]] ]}, {i,1,Length[points]}] ],
		Axes -> Automatic, opts ]
  ];

unitDirField[f_,g_,{x_,xmin_,xmax_,xnum_},
	{y_,ymin_,ymax_,ynum_}, opts___]:=
Module[{points,values,linelength,i,j,dx,dy},
	dx = N[(xmax - xmin)/(xnum-1)] ;
	dy = N[(ymax - ymin)/(ynum-1)] ;
	linelength = Min[dx/2, dy/2] ;
	points=Flatten[Table[{i,j}, 
		{i,xmin,xmax,dx},{j,ymin,ymax,dy}], 1];
	values=Flatten[Table[{f[i,j],g[i,j]}, 
		{i,xmin,xmax,dx},{j,ymin,ymax,dy}], 1];
	Show[Graphics[
		Table[Arrow[linelength,
		Theta[{0,0},values[[i]]], points[[i]] ], 
		{i,1,Length[points]}] ],
		Axes -> Automatic, opts ]
  ];
  
End[];

EndPackage[];
On[General::spell1];

Notebook Overview

A direction field is simply a prescription that tells you which direction to move given the point you're at.  A system of differential equations of the form

dx = f[x,y] dt
	dy = g[x,y] dt

Tells the velocity vector {dx/dt , dy/dt} given the values of  x  and  y  needed to compute  f[x,y]  and  g[x,y].  A plot of segments in the direction of the vector  {f[x,y] , g[x,y]}  tells the direction of motion of the solutions to the differential equations.

References to the Text

Section 22.2

In Section 22.2 we study autonomous equations in two dimensions.  This involves a growth law that does not depend on time.  It invovles the analysis of growth just by direction change.  In this analysis direction fields are used.

Definition 22.3 Direction Field

An x-y-plane with arrows pointing in the direction of x-y-change corresponding to an autonomous differential equation is called a two dimensional direction field.

Definition 22.4 Equilibrium Point

A point (xe,ye) on a direction field for the vector differential equation

dx = f(x,y) dt
	dy = g(x,y) dt

with f(xe,ye) = g(xe,ye) = 0 is called an equilibrium point.  An equilibrium point is called a local attractor if solutions that start from nearby initial conditions tend toward the point.  An equilibrium point is called a repeller if solutions that start from nearby initial conditions tend away from the point.

Sample Computations

2-D Direction Fields

The next input cell has a sample computation with our direction field function defined in the Initialization cell above.

Arrows are shaded so dark means long or strong, light means short:

Clear[x, y, f, g]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); shadeDirField[f, g, {x, -2, 10, 20}, {y, -2, 10, 20}, 
      AspectRatio -> 1];

Un-shaded arrows:

unitDirField[f, g, {x, -2, 10, 20}, {y, -2, 10, 20}, 
      AspectRatio -> 1];

Mathematica's VectorField Package

Mathematica has a vector field package with many arrow and graphics options.  You may want to investigate it further.

Needs["Graphics`PlotField`"]PlotVectorField[{f[x, y], g[x, y]}, {x, -2, 10}, {y, -2, 10}, 
      PlotPoints -> 20];

Practice Exercise

Plot the direction field associated with the differential equations

dx = x(1 - x/7 - y/7)
	dy = y(1-x/10 - y/5)

Exercises

Exercise 22.2.1 Range Wars: Cows & Sheep

Plot the direction field associated with the differential equations

dx = 2 x (1 - (x/4) - (y/4)) dtdy = y (1 - (x/3) - (y/6)) dt

Compare your solution to the CowSheep example in the Flow2D program.

Problem22.1 Predator-Prey: Rabbits & Foxes

The Lotka-Volterra equations describe a simple predator-prey system.  The equations are:

dx = (p - m y) x dt
	dy = (-q + n x) y dt

One of these variables represents the density of rabbits, and the other represents the density of foxes in an environment with unlimited space and grass.

Analyze the cases where one species is absent (x = 0 and y = 0) as well as how one species effects the other.  Which variable is foxes and which rabbits ?

Sketch the lines where the vector field is horizontal.  Which population has leveled off along these lines ?

Sketch the lines where the vector field is vertical.  Which population has leveled off along these lines ?

Make a direction field and use it to hand sketch the flow of the equation.

Check your work with this notebook when p = 0.1, q = 0.04, m = 0.005, n = 0.00004.  These are realistic values of the parameters, but make the problem badly scaled on the computer.  There have to be lots of bunnies to feed those hungry foxes, but then we have big numbers on one axis and small ones on the other.  Check the solutions:

p = 0.1; q = 0.04; m = 0.005; n = 0.00004; f[x_, y_] := (p - m y) x; g[x_, y_] := (-q + n x) y; Solve[f[x, y] == 0]Solve[g[x, y] == 0]

One Approach to the Scaling Problem

One computer solution is to replace  x  and  y  by properly scaled variables such as  u = x/1000  and  v = y/20.  Then

du = dx/1000 = (1/10 - y 5/1000)x/1000 dt
	dv = dy/20 = (-4/100 + x 4/100000)y/20 dt

du = (1/10 - v/10)u dt
	dv = (-1/25 + u/25)v dt

du = (1/10)(1 - v) u dt
	dv = (1/25)(-1 + u) v dt

Plot both the poorly scaled problem with the original  f[x,y] and g[x,y] and the re-scaled problem where  u  and  v  are fractions of the equilibrium values of rabbit and fox populations.

Compare your solution to the flow in the FoxRabbit example of the Flow2D program.

Clear[u, v, f, g]; f[u_, v_] := (1/10) u (1 - v); g[u_, v_] := (1/25) (u - 1) v; shadeDirField[f, g, {u, -1, 4, 20}, {v, -1, 4, 20}, 
      AspectRatio -> 1]; unitDirField[f, g, {u, -1, 4, 20}, {v, -1, 4, 20}, 
      AspectRatio -> 1];

Problem 22.2 Escaping Infection

The S-I-R differential equations are:

ds = - a s i  dt
	di = (a s i - b i) dt

for  s = the susceptible fraction of the population, 0 < s < 1, and  i = the infctious fraction of the population, 0 < i < 1.  Where is this vector field horizontal (in terms of  a  and  b.0)?  Is it ever vertical?

Plot the direction field in the case a = 3, b = 1.  Compare this solution to a flow animation of the solution using the Flow2D program.

Exercise 24.3.4 Classifying Equilibria

The differential equations

dx = y(1 -x/10 - y/5) dt
	dy = x(1 - x/7 - y/7) dt

have equilibria at the points (0,0), (0,5), (7,0), (4,3).  Plot its direction field to get an idea of the behavior near each of the equilibria.  Compare this solution to a flow animation of the solution using the Flow2D program.