2 Dimensional Flows
for
dx = f[x,y] dt
dy = g[x,y] dt

by
K. D.  Stroyan
University of Iowa

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

Special Mathematica Functions: flow2D[.] & vectorFlow2D[.]

Off[General::spell];
Off[General::spell1];
Off[InterpolatingFunction::dmval];
Off[ParametricPlot::pptr];
Off[Graphics::gprim];
Off[Graphics::gptn];
Off[NDSolve::ndsz];

Needs["Graphics`Colors`"]
	
colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";
Norm[a_?VectorQ] := Sqrt[a.a];

tol = .5*10^-6;

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (.075 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

		
uArrow[vec_,opts___]:= 
	Module[{a,b},
		If[Norm[N[vec]]<1,b = vec,b=vec/Norm[vec]];
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+b,opts]]
				];

flow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

vectorFlow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
								
	arr[tt_] := Map[uArrow[{f[#[[1]],#[[2]]],g[#[[1]],#[[2]]]},Tail->#,Red]&,
	{x[t],y[t]}/.soln/.t->tt];
	
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],arr[dt/100]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],arr[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

BeginPackage["Graphics`flow2D`"]
flow2D::usage =
"flow2D[{f,g},{x,y},initialValues,tFinal,dt,(options)]
	shows the flow solution of dx = f[x,y] dt, dy = g[x,y] dt 
	starting with the initial (x,y)-values = initialValues = {{a,b},{c,d},...}
	flowing up to time tFinal with animation frames every dt time units."

Begin["`Private`"]
stuff goes here 
End[]
EndPackage[]

Notebook Overview

A pair of autonomous differential equations (as above where  f[x,y]  and  g[x,y]   do not depend on time) can be viewed as a "velocity field" of "vector field."  If we think of  (x,y)  as the position in the plane and of  t  as the time, then the vector  (dx/dt , dy/dt)  is the velocity vector.  We can compute this velocity from the position  (x,y)  using the functions  (f[x,y],g[x,y]).  You might think of this as knowing the velocity vectors on the surface of a smoothly flowing river.  The differential equation tells you the velocities without telling you the motion paths.

Motion along this  "velocity field" requires that our time-path  (x[t],y[t])  satisfies  (x'[t],y'[t]) = (f[x[t],y[t]],g[x[t],y[t]]).  You can think of a solution as the actual motion of a particle on the surface of the river.  Each particle begins at a different point and the flow of the surface is the aggregate motion of all the solutions with different initial conditions.

This NoteBook contains a package that will trace the motion of any number of particles.  You tell it the "velocity field"  (f[x,y],g[x,y])  and the intial positions of the particles you wish to follow.  The package builds an animation of the motion.

References to the Text

This NoteBook goes with Chapter 22 of the text.

Step-by-step Computation

This section gives a step-by-step explanation of the computations of this NoteBook.  Once you understand the steps, you may want to combine them in a single longer computation.

The Flows of Text Example 22.3: Competition

To use the flow2D[] package which you <Enter>ed above as the initialization, you need to:

(1)  Define two functions, f[x,y] and g[x,y], as below that do not depend on time.
(This makes the equations autonomous, so the flow idea applies.)

(2)	Enter a list of initial {x,y} pairs, {{a,b},{c,d},{e,f},..}.

(3)	Give the final time and time step between animation frames.

(4)	Enter this information in the flow2D[] function
The following does each of these steps in a separate cell.

The functions of the differential equation

Clear[f, g, x, y, t]; f[x_, y_] := x (1 - (x/7) - (y/7)); g[x_, y_] := y (1 - (x/10) - (y/5)); f[x, y]g[x, y]

The initial conditions

initVals = {{10, 1}, {10, 5}, {10, 10}, {1, 9}, {5, 9}, {1, 1}}

The final time and time step between graphs

tFinal = 10; dt = 1; Print[Floor[(tFinal/dt)], " graphs from t = 0 to t = ", tFinal]

The flow

To show the flow as a real time animation, highlight the cell containing the graphs and press the <Command> + <y>  keys at the same time. "vectorFlow2D[.]" shows the point moving under the flow produced by the vector field and shows the direction vector as you go.  "flow2D[.]" only shows the point moving.

solutions = 
      vectorFlow2D[{f, g}, {x, y, t}, initVals, tFinal, dt, 
        AxesLabel -> {"x", "y"}, 
        PlotRange -> {{-1, 10}, {-1, 10}}];

View the Solutions

ParametricPlot[Evaluate[{x[t], y[t]} /. solutions], {t, 0, tFinal}, 
      AspectRatio -> 1, PlotRange -> {{-1, 10}, {-1, 10}}, 
      AxesLabel -> {"x", "y"}]; Plot[Evaluate[x[t] /. solutions], {t, 0, tFinal}, 
      PlotRange -> {0, 10}, AxesLabel -> {"t", "x"}]; 
    Plot[Evaluate[y[t] /. solutions], {t, 0, tFinal}, 
      PlotRange -> {0, 10}, AxesLabel -> {"t", "y"}]; 
    Print["Solutions at time 5"]; Evaluate[{x[t], y[t]} /. solutions /. t -> 5.0]Print["Solutions at time 10"]; Evaluate[{x[t], y[t]} /. solutions /. t -> 10.0]

Combined Computation

This section gives the complete computation of this NoteBook in a single computation.  This is more convenient than <Entering> the different portions of the computation in separate cells.

The Flow of Text Example 22.4: Fierce Competition

Be patient, this computation takes a long time, because there are a lot of initial conditions and a lot of time steps.

Clear[f, g, x, y, t]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); initconds = 
      Flatten[Table[{x, y}, {x, - .1, 9.9, 2}, {y, - .1, 9.9, 2}], 
        1]; tfinal = 10; dt = 1; Print["Flow solution of the differential equations:"]; Print["dx = ", f[x, y], " dt"]; Print["dy = ", g[x, y], " dt"]; flow2D[{f, g}, {x, y, t}, initconds, tfinal, dt, 
      AxesLabel -> {"x", "y"}, 
      PlotRange -> {{-1, 10}, {-1, 10}}];

To show the flow as a real time animation, highlight the cell containing the graphs and press the <Command> + <y>  keys at the same time.

Exercises

Exercise 22.2.1 (Range wars: Cows and sheep compete for grass)

Modify the change functions to see the flow of the equations:

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

All you really need to do is re-define the functions  f[x,y]  and  g[x,y].

\(Clear[f, g, x, y, t]; 


    f[x_, y_]\  := \ 
      \(?? ;
g[x_, y_]\  := \ 
          \(?? ;


                initxys = {{1, 1}, {1, 2}, {2, 1}, {3, 3}, {1, 3}, {3, 
                  1}, { .5, 1}, {1,  .5}, {2.5, 3}, {3, 2.5}, {3.5, 2.5}}; 

            \(tfinal = 10; \)
\(dt = 1; \)

              \(Print["Flow solution of the competing species differential equations:"]; \)
\(Print["dx = ", f[x, y], " dt"]; \)

              \(Print["dy = ", g[x, y], " dt"]; \)

              \(vectorFlow2D[{f, g}, {x, y, t}, initxys, tfinal, dt, 
                PlotRange -> {{\(- .1\), 6}, {\(- .1\), 6}}, 
                AxesLabel -> {cows, sheep}]; \)\)\)\)

Predator-prey interactions: Foxes Eat Bunnies

Problem 22.1

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 compass heading chart and a direction field and use them to hand sketch some parametric solutions of the equation.

Check your work with the FoxRabbit example in Flow2D when p = 0.1, q = 0.04, m = 0.005, n = 0.00004.

After watching the flow animation you should sketch the solutions as functions of time, x vs. t and y vs. t.

Use this notebook to solve the equations with the initial conditions (1000,5), (1000,20). Compare these solutions.

Copy computations above and modify the functions to see the flow of the equations.   In this example, you will probably want to use the initial conditions:  initconds = Table[{1000,y},{y,5,45,5}];
(This problem is poorly scaled, so we drop the arrows on the flow.  See the DirField predator-prey example for a re-scaling.)

Clear[f, g, x, y, t];

f[x_, y_] :=  ??;
g[x_, y_] :=  ??;

\(p = 0.1; 
\(q = 0.04; 
              \)\)
m = 0.005; 
n = 0.00004; 
initxys = Table[{1000, y}, {y, 5, 45, 5}]; 
tfinal = 100; 
dt = 10; 
Print["Flow solution of the predator-prey differential \
equations:"]; 
Print["dx = ", f[x, y], " dt"]; 
Print["dy = ", g[x, y], " dt"]; 
flow2D[{f, g}, {x, y, t}, initxys, tfinal, dt, 
                PlotRange -> {{-1, 4400}, {-1, 44}}, 
                AxesLabel -> {Rabbits, Foxes}];

Exercise CD 24.4.2: The Lotka-Volterra Invariant

invar := p Log[y] + q Log[x] - m y - n x; ContourPlot[invar, {x, 10, 4000}, {y, 0.01, 40}];

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 flow in the case a = 3, b = 1.

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 flow animation using the Flow2D program.

Check the four local results by having enough initial conditions near each equilibrium.