3-D Flow of the Autonomous System
dx = f[x,y,z] dt
dy = g[x,y,z] dt
dz = h[x,y,z] 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]]
				];
				
flow3D[{f_,g_,h_},{x_,y_,z_,t_},initVals_,tFinal_,dt_,opns___] :=
Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
NDSolve[{x[0]==#[[1]],y[0]==#[[2]],z[0]==#[[3]],
	x'[t] == f[x[t],y[t],z[t]],
	y'[t] == g[x[t],y[t],z[t]],
	z'[t] == h[x[t],y[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},opns,PlotRange->{{0,10},{0,10},{0,10}},
DisplayFunction->Identity];
			
pt[tt_] := Map[Graphics3D[{Blue,AbsolutePointSize[4],
Point[{#[[1]],#[[2]],#[[3]]}]}]&,
{x[t],y[t],z[t]}/.soln/.t->tt];

Show[{crv[dt/100],pt[0]},opns,
PlotRange->{{0,10},{0,10},{0,10}},AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
Return[soln];
];

vectorFlow3D[{f_,g_,h_},{x_,y_,z_,t_},initVals_,tFinal_,dt_,opns___] :=
Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
NDSolve[{x[0]==#[[1]],y[0]==#[[2]],z[0]==#[[3]],
	x'[t] == f[x[t],y[t],z[t]],
	y'[t] == g[x[t],y[t],z[t]],
	z'[t] == h[x[t],y[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},
opns,PlotRange->{{0,10},{0,10},{0,10}},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]],#[[2]],#[[3]]],
g[#[[1]],#[[2]],#[[3]]],h[#[[1]],#[[2]],#[[3]]]},
Tail->#,Red]&,{x[t],y[t],z[t]}/.soln/.t->tt];

Show[{crv[dt/100],pt[0],arr[dt/100]},opns,
PlotRange->{{0,10},{0,10},{0,10}},AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt],arr[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 once colors woik
End[]
EndPackage[]

Notebook Overview

This notebook will perform three dimensional flows of an autonomous system.

References to the Text

This NoteBook goes beyond Chapter 22 of the text.

Four Views of the Explicit Solution of a 2-D System

The following computations show different views of the plot of the vector function  X[t] = {x[t],y[t]} vs the independent vaariable t.  This is done with the 3-D Flow of
	
		dx = f[x,z] dt
		dy = dt
		dz = g[x,z] dt

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

f[x_,y_,z_] := x(1 - x/7 - z/7);
h[x_,y_,z_] := z(1 - x/10 - z/5);
g[x_,y_,z_] := 1;

initVals = {
		{10,0,1},{10,0,5},{10,0,10},
		{1,0,9},{5,0,9},
		{1,0,1}
		};
tFinal = 10 ;
dt = 1 ;
Print["A balanced view of solution with time:"];
solutions=vectorFlow3D[ {f,g,h}, {x,y,z,t}, initVals,
		tFinal, dt, ViewPoint->{5,3,4},AxesLabel->{"x","t","y"},
		Ticks->None,
		PlotRange -> {{0,10},{0,10},{0,10}} ];

Print["Seeing mostly y = y[t] vs x = x[t]"];
vectorFlow3D[ {f,g,h}, {x,y,z,t}, initVals,
		tFinal, dt, ViewPoint->{-1,18,1},AxesLabel->{"x","t","y"},
		Ticks->None,
		PlotRange -> {{0,10},{0,10},{0,10}} ];

Print["Seeing mostly y = y[t] vs t"];
vectorFlow3D[ {f,g,h}, {x,y,z,t}, initVals,
		tFinal, dt, ViewPoint->{18,-1,-1},AxesLabel->{"x","t","y"},
		Ticks->None,
		PlotRange -> {{0,10},{0,10},{0,10}} ];

Print["Seeing mostly x = x[t] vs t"];
vectorFlow3D[ {f,g,h}, {x,y,z,t}, initVals,
		tFinal, dt, ViewPoint->{1,-1,-18},AxesLabel->{"x","t","y"},
		Ticks->None,
			PlotRange -> {{0,10},{0,10},{0,10}} ];

ParametricPlot[Evaluate[{x[t],z[t]}/.solutions],
{t,0,tFinal},AspectRatio->1];
Plot[Evaluate[x[t]/.solutions],{t,0,tFinal}];
Plot[Evaluate[z[t]/.solutions],{t,0,tFinal}];
Plot[Evaluate[y[t]/.solutions],{t,0,tFinal}];

Practice Exercises

Combined Example

The commands in the cell below will perform a three dimensional flow of an autonomous system.

Clear[f,g,h,x,y,z,t];
a = 2;
a1 = .5;
a2 = .5;
a3 = .5;
b = 1; c = 1; e = 1;
f[x_,y_,z_] := a y - a1 x - a2 x ;
g[x_,y_,z_] := a1 x - a3 y - b z y ;
h[x_,y_,z_] := -c z + e z y ;

initconds = Flatten[ Table[ {x, y, z}, {x, 0, 5}, {y, 0, 5},{z, 0, 5} ], 2 ];

flow3D[{f, g, h},{x,y,z,t}, initconds, 1, .2,
	PlotRange->{{0,10},{0,10},{0,10}},AxesLabel->{"x","y","z"}];

Three Competing Species

Clear[f,g,h,x,y,z,t];
f[x_,y_,z_] := x(1-x/5-y/12-z/21);
g[x_,y_,z_] := y(1-x/15-y/6-z/14);
h[x_,y_,z_] := z(1-x/10-y/18-z/7);
initconds = Flatten[
		Table[{a,b,c},{a,1,9,4},{b,1,9,4},{c,1,9,4}],2];
flow3D[{f, g, h},{x,y,z,t},initconds, 20, 2,
	PlotRange->{{0,10},{0,10},{0,10}},AxesLabel->{"x","y","z"}];

Exercises

None.  The text does not have 3-D dynamical systems.