Various Views of Solutions
of
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 used in this NoteBook

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;
								
OrthUp[a_?VectorQ] :=
	Module[{b,c,d,f},
		d = N[a/Norm[a]];
		If[Chop[d[[3]]^2 - 1,tol] == 0, b={1,0,0} ,
		b = {0,0,1} - d[[3]] d; b = b/Norm[b];];
	Return[b]];
	
OrthSet[a_?VectorQ] :=
	Module[{b,c,d},
		d = N[a/Norm[a]];
		If[Chop[d[[1]],tol]==0,b={1,0,0},
			 If[Chop[d[[2]],tol]==0,b={0,1,0},
			  b={-d[[2]],d[[1]],0}]];
			b=b/Norm[b];c = Cross[d,b];
	Return[{b,c/Norm[c],d}]];

				
Rotate[a_?VectorQ,theta_] :=
	Module[{u,v,r},
		u = OrthSet[a];
		v = Transpose[u];
		r = N[{{Cos[theta],-Sin[theta],0},{Sin[theta],Cos[theta],0},{0,0,1}}];
	Return[v.r.u]];


Vector3D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,u2,u3,u4,rotmat,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[Graphics3D[Point[{0,0,0}]]];					
	,		
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics3D[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (ntip-ntail);
	rotmat = Rotate[Punit,N[Pi/2]];
	up = OrthUp[Punit];
	Punit = .075 Punit;
	up = .15 Norm[Punit] up;
	u2 = rotmat.up;
	u3 = rotmat.u2;
	u4 = rotmat.u3;
	arrow = Line[{ntail,ntip,
				ntip - Punit + up,ntip,
				ntip - Punit + u2,ntip,
				ntip - Punit + u3,ntip,
				ntip - Punit + u4}];
	]];
Return[Graphics3D[{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,0}};
				Return[Vector3D[a,a+b,opts]]
				];


solnViews[{f_,h_},{x_,z_,t_},initVals_,tFinal_,dt_,opns___] :=
Module[{soln,arr,pt,crv,tt,y,flo3,xPlot,yPlot,phasePlot},
soln = Flatten[Map[
NDSolve[{x[0]==#[[1]],y[0]==0,z[0]==#[[2]],
	x'[t] == f[x[t],z[t]],
	y'[t] == 1,
	z'[t] == h[x[t],z[t]]},{x[t],y[t],z[t]},{t,0,tFinal}]&,
initVals],1];
										
crv[tt_] :=
ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.soln],{t,0,tt},
PlotRange->{{0,10},{0,10},{0,10}},Ticks->None,
DisplayFunction->Identity];
			
pt[tt_] := Map[Graphics3D[{Blue,AbsolutePointSize[4],
Point[{#[[1]],#[[2]],#[[3]]}]}]&,
{x[t],y[t],z[t]}/.soln/.t->tt];

arr[tt_] := Map[uArrow[{f[#[[1]],#[[3]]],
1,h[#[[1]],#[[3]]]},
Tail->#,Red]&,{x[t],y[t],z[t]}/.soln/.t->tt];

flo3[tt_] := Show[{crv[tt],pt[tt],arr[tt]}, ViewPoint->{5,3,4},
DisplayFunction->Identity];

xPlot[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t, 0, tt},opns,
AspectRatio->1,AxesLabel->{"x","t"},DisplayFunction->Identity];

yPlot[tt_] := ParametricPlot[Evaluate[{y[t],z[t]}/.soln],{t, 0, tt},opns,
AspectRatio->1,AxesLabel->{"t","y"},DisplayFunction->Identity];

phasePlot[tt_] := Show[
ParametricPlot[Evaluate[{x[t],z[t]}/.soln],{t, 0, tt},opns,
AspectRatio->1,AxesLabel->{"x","y"},DisplayFunction->Identity],
Graphics[{Blue,AbsolutePointSize[4],Point[Flatten[Evaluate[{x[t],z[t]}/.soln/.t->tt],1]]}]];

Do[
Show[GraphicsArray[{{xPlot[tt],flo3[tt]},{phasePlot[tt],yPlot[tt]}}],
DisplayFunction->$DisplayFunction];
,{tt,dt,tFinal,dt}];
Return[soln];
];

Notebook Overview

This NoteBook shows both the explicit plots x[t] vs. t and y[t] vs. t and the parametric plot {x[t],y[t]} of the solution of an initial value problem

x[0] = x0
	y[0] = y0
	dx = f[x,y] dt
	dy = g[x,y] dt

The "flow" of the differential equation is an animated parametric plot {x[t],y[t]} and the geometric analysis of the solutions of the initial value problems associated with the differential equations (or vector field {f[x,y],g[x,y}}) rely on this parametric representation.  It is important for you to see the conncetion between the more familiar explicit plots and the parametric plot.

References to the Text

This NoteBook goes with Chapter 22 of the text.

Step-by-step

We begin with another look at the dynamical system of the fierce competition example 22.?? of the text and the combined computation of the Flow2D NoteBook. <Enter> the next computation to see the dead plot of the solutions that were animated in Flow2D.

Clear[f, g, x, y, t]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); initXYs = Flatten[Table[{x, y}, {x, 0, 10, 2}, {y, 0, 10, 2}], 1]; tFinal = 10; soln=Flatten[(NDSolve[{x[0] == #1[[1
                              ]],y[0] == #1[[2
                              ]],x'[t]==f[x[t], y[t]],y'[t]==g[x[t], y[t]]},{x[t], y[t]},{t, 0, tFinal}]&)/@initXYs,1];ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, tFinal}, 
      AspectRatio -> 1, AxesLabel -> {"x", "y"}];

This dynamical system has 3 equilibria.  One at (0,7) on the y-axis, one at (10,0) on the x-axis, and a third at (4,3).  Solutions typically begin heading toward the central equilibrium and then veer off and head toward one to the two axis attractors. (Unless they hit (4,3) exactly.)

The solution starting at (2,2)

Consider the solution to
	x[0] = 2,
	y[0] = 2,
	dx = x(1 - x/10 - y/5) dt,
	dy = y(1 - x/7 - y/7) dt
by <Enter>ing the next computation.  This will show the parametric plot starting at {2,2}, the plain x[t] plot and the plain y[t] plot:

Clear[f, g, x, y, t]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); tFinal = 10; soln=NDSolve[{x[0] == 2,y[0] == 2,x'[t]==f[x[t], y[t]],y'[t]==g[x[t], y[t]]},{x[t], y[t]},{t, 0, tFinal}];ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, tFinal}, 
      AspectRatio -> 1, AxesLabel -> {"x", "y"}, 
      PlotRange -> {{0, 10}, {0, 10}}]; Plot[Evaluate[x[t] /. soln], {t, 0, tFinal}, AxesOrigin -> {0, 0}, 
      PlotRange -> {0, 4}, AxesLabel -> {"t", "x"}]; Plot[Evaluate[y[t] /. soln], {t, 0, tFinal}, AxesOrigin -> {0, 0}, 
      PlotRange -> {0, 7}, AxesLabel -> {"t", "y"}];

These plots show x first increasing towart the central equilibrium, then decreasing back toward zero and y steadily increasing.  However, the plots are annoying because of the different scales.  The next approach puts things on a more equal scaling.

A Combined Look at the Solution

In order to see the connection between the different plots on the same scales and in a more convenient orientation, we "solve fot t" by considering it the solution to the initial value problem

z[0] = 0	
	dz = dt

and adding that to our original problem.  <Enter> the next cell now to see how this looks:

Clear[f, g, x, y, z, t]; f[x_, y_] := x (1 - (x/10) - (y/5)); g[x_, y_] := y (1 - (x/7) - (y/7)); tFinal = 10; soln=NDSolve[{x[0] == 2,y[0] == 2,z[0] == 0,x'[t]==f[x[t], y[t]],y'[t]==g[x[t], y[t]],z'[t]==1},{x[t], y[t], z[t]},{t, 0, tFinal}];xtPlot[tt_] := 
      ParametricPlot[Evaluate[{x[t], z[t]} /. soln], {t, 0, tt}, 
        AspectRatio -> 1, AxesLabel -> {"x", "t"}, 
        PlotRange -> {{0, 10}, {0, 10}}, 
        DisplayFunction -> Identity]; txPlot[tt_] := 
      ParametricPlot[Evaluate[{z[t], x[t]} /. soln], {t, 0, tt}, 
        AspectRatio -> 1, AxesLabel -> {"t", "x"}, 
        PlotRange -> {{0, 10}, {0, 10}}, 
        DisplayFunction -> Identity]; yPlot[tt_] := 
      ParametricPlot[Evaluate[{z[t], y[t]} /. soln], {t, 0, tt}, 
        AspectRatio -> 1, AxesLabel -> {"t", "y"}, 
        PlotRange -> {{0, 10}, {0, 10}}, 
        DisplayFunction -> Identity]; phasePlot[tt_] := 
      ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, tt}, 
        AspectRatio -> 1, AxesLabel -> {"x", "y"}, 
        PlotRange -> {{0, 10}, {0, 10}}, 
        DisplayFunction -> Identity]; Show[GraphicsArray[{{xtPlot[tFinal], txPlot[tFinal]}, {phasePlot[tFinal], 
            yPlot[tFinal]}}], DisplayFunction -> $DisplayFunction];

Notice that the top row of this graphic has two plots of  x[t] vs. t.  The conventional one at the right and a "backwards one at the left directly over the parametric plot.  If you follow a point along the parametric plot, the point on the "backwards" explicit x[t] vs t graph lies directly above this point.

Similarly, if you follow a point along the parametric graph, the y[t] vs. t point appears on the graph to the right in the bottom row.

Practice Exercise

Plot the solution to these same differential equations starting at th point {4,8} and describe in words how the x[t] and y[t] solutions increase or decrease as a point moves along the parametric plot.

Combined Computation

We could graph the vector {x[t],y[t]} explicitly using three dimensions.  In this case, the parametric plot {x[t],y[t]} is just a "top" view of the solution, the x[t] vs. t plot an "end" view and the y[t] vs. t solution a "side" view.  This can all be easily animated as shown next for a well-scaled predator-prey system

x[0] = 1/2,
	y[0] = 1.0,
	dx(hares) = x(3-y),
	dy(foxes) = (x-2)y

Clear[hares, foxes, x, y, t]; hares[x_, y_] := x (3 - y); foxes[x_, y_] := (x - 2) y; initVals = {{ .5, 1}}; tFinal = 10; dt =  .5; solutions = 
      solnViews[{hares, foxes}, {x, y, t}, initVals, tFinal, dt, 
        PlotRange -> {{-1, 10}, {-1, 10}}];

Practice Exercise

Animate the {hares,foxes} dynamics when x[0] = 2 and y[0] = 2.

Exercises

Exercise 22.3.1

Review the 2-D flow of a basic system such as the Cow-Sheep-equations

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

Then make an animation of the various views of their solutions first with the initial conditions

x[0] = 2, y[0] = 1
then
	x[0] = 1, y[0] = 2

and some other values you find interesting.(Can you find initial conditions that hit the central equilibrium?)

Describe the way x[t] and y[t] increase or decrease as the parametric plot moves toward equilibrium.