Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 30: Continuous Production and Exchange

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

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[]

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Get Started

Read the printed project and solve Hints 30.1, 2, 3.  Following is help with computing.

Computer Experiments

Define the basic differential equations:

\(Clear[f, g, x, y, t, p, q, r, c, λ, Subscript]; 

    f[x_, y_] := c\ \((p\/\(p\ x + q\ y + 1\) - r)\); 

    g[x_, y_] := c\ \((p\/\(q\ x + p\ y + 1\) - r)\); 

    Print["Flow solution of the differential equations:"]; 

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

Find the general equilibrium point:

\(soln = Solve[{f[x, y] == 0, g[x, y] == 0}, {x, y}]; 
    
{x\_e, y\_e} = {x, y} /. soln[\([1]\)]\)

Find the characteristic polynomial coefficients as int hte text Procedure 24.5

aX = D[f[x, y], x] /. {x -> (Subscript[x, e]), y -> (Subscript[y, e])}aY = D[f[x, y], y] /. {x -> (Subscript[x, e]), y -> (Subscript[y, e])}bX = D[g[x, y], x] /. {x -> (Subscript[x, e]), y -> (Subscript[y, e])}bY = D[g[x, y], y] /. {x -> (Subscript[x, e]), y -> (Subscript[y, e])}\(soln = 
      Solve[λ^2 - \((aX + bY)\)\ λ\  + 
            \((aX\ bY - aY\ bX)\) == 0, λ]; 

    λ\_1 = λ /. soln[[1]]\)(Subscript[λ, 2]) = λ /. soln[[2]]

Then look at the following for various cases.

First Test: p  = 7/8, r = 1/7

\(c\  = \ 1; 
q = 1 - p; 
p\  = \ 7/8; 
r = 1/7; 

    initconds = 
      Union[Table[{x, 1}, {x, 1, 9, 2}], Table[{1, y}, {y, 1, 9, 2}]]; 

    initconds = Union[initconds, Table[{9, y}, {y, 1, 9, 2}]]; 

    initconds = Union[initconds, Table[{x, 9}, {x, 1, 9, 2}]]; 

    tfinal = 100; 
dt = 25; 

    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}}]; 

    Print[\*""Equilibrium at {(Subscript[x, e]),(Subscript[y, e])} ="", 
      N[{x\_e, y\_e}]]; 

    Print[\*""Characteristic roots from text 24.5, {(Subscript[λ, 1]),(Subscript[λ, 2])}="", 
      N[{λ\_1, λ\_2}]]; \)

Conclusions about this example?

First Test: p  = 1/3, r = 1/7

\(c\  = \ 1; 
q = 1 - p; 
p\  = \ 1/3; 
r = 1/7; 

    initconds = 
      Union[Table[{x, 1}, {x, 1, 9, 2}], Table[{1, y}, {y, 1, 9, 2}]]; 

    initconds = Union[initconds, Table[{9, y}, {y, 1, 9, 2}]]; 

    initconds = Union[initconds, Table[{x, 9}, {x, 1, 9, 2}]]; 

    tfinal = 100; 
dt = 25; 

    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}}]; 

    Print[\*""Equilibrium at {(Subscript[x, e]),(Subscript[y, e])} ="", 
      N[{x\_e, y\_e}]]; 

    Print[\*""Characteristic roots from text 24.5, {(Subscript[λ, 1]),(Subscript[λ, 2])}="", 
      N[{λ\_1, λ\_2}]]; \)

Conclusions about this example?

Conclusions on the Experiments

Finish The Theory with Hints 30.6, 7

Read these hints and solve the problem.  Compare with your experiments.