DE.04 Modern DiffEq Issues
TUTORIALS

Needs["Graphics`FilledPlot`"];

If[MemberQ[{"Global`","Graphics`Common`GraphicsCommon`","Vector3D`"},Context[ScaleFactor]],Remove["*`ScaleFactor"]];
<<"Graphics`Legend`";
<<"Graphics`Common`GraphicsCommon`";
Remove["*`ScaleFactor"];

If[MemberQ[{"Global`", "Miscellaneous`SIUnits`", "Graphics`Colors`"},Context[Gray]], Remove["*`Gray"]];
<<"Miscellaneous`Units`";
Remove["*`Gray"];
<<"Graphics`Colors`";

Graphics`Colors`GosiaGreen=RGBColor[0, 0.392187, 0];
Off[General::spell];
Off[General::spell1];
Off[ParametricPlot3D::ppcom];
Off[ParametricPlot::ppcom];
Off[Plot::plnr];
Needs["Calculus`Master`"];
Off[Unique::"usym"];
SetOptions[Limit,Analytic->True];

CMView = {2.7, 1.6, 1.2};

Off[Graphics::gpt];

New Vector Drawers

(* :Title:  NewVectorDraw  *)
(* :Author: Bruce Carpenter *)
(* :Summary:  Construct 2 and 3 dimensional vectors,
     arrows, and arrowheads as graphics 
     objects.
*)
(* :Date:  Copyright 1997, Math Everywhere, Inc. 
*)
(*  :Drawbacks:  Since the vectors are drawn out of 
     context, the shapes of the arrowheads are only 
     correct with AspectRatio set to Automatic.
*)
BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2
Perpend[a_?VectorQ] := Normalize[If[a[[1]]==0,{1,0,0},{0,a[[3]],-a[[2]]}]]/;Length[a]==3

Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

VectorField and Phase Lines

(*  :Title: VectorAndPhase.m  *)

(*  :Author:  Bill Davis  *)

(*  :Dates:  about March, 1996 *)

(*  :Drawbacks:  I wish someone serious would come and turn these into
                 Really carefully done packages.  *)

BeginPackage["PhaseStuff`","Vector3D`",
		"Utilities`FilterOptions`"];

PhaseLine::usage = "PhaseLine[fcn,{var,a,b},tint,pos]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint,  at the x-coordinate pos.";
	
PhaseLine3D::usage = "PhaseLIne3D[fcn,{var,a,b},tint,loc]
	makes a vertical phase line for the function, fcn, of the variable, var,
	for the variable going from a to b, with 
	a color given by tint, and base point at loc = {x,y}";
	
VectorField::usage = "VectorField[{function,{var1,a,b,space1},{var2,c,d,space2}]
	will make a vector field of normalized vectors, colored by
	magnitude, for the function of the variables on the
	range {r1,r2} with spacing between points given by 
	{space1, space2}.  If a single function of a single variable is given,
	the vector field {1,f} will be plotted.";
	
StikPikXZ::usage = "StikPikXZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikYZ::usage = "StikPikYZ[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
StikPikXY::usage = "StikPikXY[picture,a:_0] takes a 2 dimensional
	Mathematica picture and embeds it on the XZ plane at y = a.  
	Some optional characteristics of the graphs are lost in the 
	process.  StikPikYZ and StikPikXZ are similar.";
	
Begin["`Private`"];	

PhaseLine[fcn_,{var_,a_,b_},tint_,pos_:0] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos,a},{pos,b}}]};
	dots = Point[{pos,#}]&/@zero;
	dots = Join[{PointSize[0.05],tint},dots];
	heads = Table[MakeHead[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),3/4 slen];
	Return[ArrowHead[{c,a(.5-s/5)+b(.5+s/5)},
		{0,s u},VectorColor->tint]]]
		
PhaseLine3D[fcn_,{var_,a_,b_},tint_,pos_:{0,0}] := 
	Module[{len,len1,k,zset,zero,zero1,signs,stick,dots,pointers},
	zset = Solve[fcn==0,var];
	zero = If[Length[zset]==0,{},var/.zset];
	zero = Sort[zero];
	zero = Select[zero,((# >= a)&&(# <= b))&];
	zero1 = Sort[Union[zero,{a,b}]];
	len = Length[zero];
	len1 = Length[zero1];
	signs = Table[Sign[fcn/.var->((zero1[[k+1]] + zero1[[k]])/2)],
		{k,1,len1-1}];
	stick = {Thickness[0.007],Line[{{pos[[1]],pos[[2]],a},
		{pos[[1]],pos[[2]],b}}]};
	dots = Point[{pos[[1]],pos[[2]],#}]&/@zero;
	dots = Join[{PointSize[0.03],tint},dots];
	heads = Table[MakeHead3D[zero1[[k]],zero1[[k+1]],pos,
		signs[[k]],b-a,tint],{k,1,len1-1}];
	Return[{Graphics3D[Join[{tint},stick,dots]],heads}]
		]
		
MakeHead3D[a_,b_,c_,s_,slen_,tint_] := Module[{u,v,w},
	u = Min[2(b-a),1/2 slen];
	Return[ArrowHead[{c[[1]],c[[2]],a(.5-s/5)+b(.5+s/5)},
		{0,0,s u},VectorColor->tint]]]
		
		
VectorField[function_List,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] := 
	Module[{i,j,U,V,UU,VV,ulow,uhigh,vlow,vhigh,maxval,minval,
		lengths,directions,locations,mlen},
	mlen = Min[sp1,sp2];
	normf[uu_,vv_] = Sqrt[function[[1]]^2+function[[2]]^2]/.
		{var1->uu,var2->vv};
	lengths = Table[N[normf[U,V]],{U,a,b,sp1},{V,c,d,sp2}];
	maxval = Max[lengths]; minval = Min[lengths];
	If[minval == maxval,maxval += 1];
	locations = Table[{uu,vv},{uu,a,b,sp1},
		{vv,c,d,sp2}];
	directions = Table[If[N[normf[U,V]]==0,{0,0},
		mlen (function/.{var1->U,var2->V})/(1.3 N[normf[U,V]])],
			{U,a,b,sp1},{V,c,d,sp2}];
	{UU,VV} = Dimensions[directions][[{1,2}]];
	Return[Table[Arrow[directions[[i,j]],
			Tail->locations[[i,j]],Aperture->1,
			HeadRatio->.4,ShaftWidth->0.008, VectorColor->
			RGBColor[(lengths[[i,j]] - minval)/(maxval-minval),0,
				(maxval -lengths[[i,j]])/(maxval-minval)]],{i,UU},{j,VV}]]
		]
		
VectorField[function_,{var1_,a_,b_,sp1_:1},{var2_,c_,d_,sp2_:1}] :=
	VectorField[{1,function},{var1,a,b,sp1},{var2,c,d,sp2}]
	
	
Needs["Utilities`FilterOptions`"]

FixOptions[a_] := Module[{picture = a, opt3},
If[Options[picture] != {}, 
	opt3 = FilterOptions[Graphics3D,#]&/@Options[picture];
	(* Fix things that might carry 2D characteristics. *)
	opt3 = opt3//. {(AxesLabel -> _) ->(AxesLabel -> None),
			(PlotRange->_) -> (PlotRange->Automatic),
			(PlotRegion -> _) -> (PlotRegion -> Automatic),
			(Ticks -> _) -> (Ticks -> Automatic)};
	picture = Apply[List,picture];
	picture[[-1]] = opt3;
	picture = Apply[Graphics3D,picture]];
	Return[picture]
	]

StikPikXZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,a,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikYZ[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{a,u,v}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]


StikPikXY[picture_,a_:0] := Module[{PIK3 = picture},
	PIK3 = PIK3//.{Graphics->Graphics3D,
		{u_?NumberQ,v_?NumberQ}:>{u,v,a}};
	PIK3 = FixOptions[PIK3];
	Return[PIK3]
		]
		
End[];

EndPackage[];
Null

ThreeAxes[u,v]

ThreeAxes[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,0},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
ThreeAxes[u_] := ThreeAxes[u,u/8];
ThreeAxes::usage = "ThreeAxes[a,b] makes a standard cartesian axis graphics object
	with x, y, and z running from -a to a, and with axis
	labels b units beyond the tips of the axes.  ThreeAxes[a]
	is ThreeAxes[a,a/8].";

Crisscross

Clear[x];
crisscross = Plot[{Cos[x],Sin[x]},{x,Pi/8,3 Pi/8},
    PlotStyle->{{RGBColor[1,0,0],Thickness[0.02]}},
    Axes->None,AspectRatio->1/3,
    DisplayFunction->Identity];

T.1) Sensitive dependence on starter data

T.1.a.i) Flow plots

Here's a simple differential equation:

\(Clear[diffeq, y, t, f, starter]; 
b\  = \ 10; 

    f[t_, y_]\  = \ y\  - \ 0.8\ t; 

    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

And a flow plot:

endtime = 6; scalefactor = 0.4; {ylow, yhigh} = {-4, 10}; \({tlow, thigh} = {0, 10}; 


    flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.8], {t, tlow, thigh, \(thigh - tlow\)\/10}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/8}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"t", "y"}, 
      AspectRatio -> 1\/2, 
      PlotLabel -> setup[[1, 1]]]; 
    \)

And see the plots of two solutions:

\(Clear[y1, y2, y3, t]; 
{starter1, starter2}\  = \ {0.79, 0.81}; 

    endtime\  = \ 7; 


    y1[t_]\  = \ 
      y[t] /. \(NDSolve[\((diffeq /. starter -> starter1)\), 
            y[t], {t, 0, endtime}]\)[\([1]\)]; 

    y2[t_]\  = \ 
      y[t] /. \(NDSolve[\((diffeq /. starter -> starter2)\), 
            y[t], {t, 0, endtime}]\)[\([1]\)]; 


    solutionplots\  = \ 
      Plot[{y1[t], y2[t]}, {t, 0, endtime}, 
\ \ \ \ \ \ 
        PlotStyle -> 
	{{Thickness[0.01], Red}}, 
        AspectRatio -> 1\/GoldenRatio, 
	\ \ PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, DisplayFunction -> Identity]; 


    Show[solutionplots, flowplot, DisplayFunction -> $DisplayFunction]; \)

What do folks mean when they talk about "sensitive dependence on starter data?"

You're looking at it.  
Take another look:

Show[solutionplots, flowplot, 
      DisplayFunction -> $DisplayFunction];

What you see are two solutions of:

diffeq

The starting points of the two solutions are almost the same; yet the two solutions are totally different in character. 
You gotta agree that small changes in the starting conditions can result in big changes in solutions. This is exactly what folks mean when they talk about senistive dependence on starting data.

T.1.a.ii)

Do all differential equations exhibit sensitive dependence on starter data?

Not all of them.
Here's one that doesn't exhibit overly sensitive dependence on starter data:

\(Clear[f, t, y]; 
b\  = \ 10; 
f[t_, y_]\  = \ \ \(-y\)\ ; 

    setup\  = \ 
      \((diffeq\  = \ {\(y'\)[t]\  == \ f[t, y[t]], y[0]\  == \ starter})
          \) // ColumnForm\)

And a flow plot:

endtime = 6; scalefactor = 0.4; {ylow, yhigh} = {-4, 4}; \({tlow, thigh} = {0, 8}; 


    flowplot = 
      Table[Arrow[{1, f[t, y]}, Tail -> {t, y}, 
          VectorColor -> Blue, ScaleFactor -> scalefactor, 
          HeadSize -> 0.6], {t, tlow, thigh, \(thigh - tlow\)\/8}, {y, 
          ylow, yhigh, \(yhigh - ylow\)\/6}]; 


    Show[flowplot, Axes -> True, AxesLabel -> {"t", "y"}, 
      AspectRatio -> 1\/GoldenRatio, 
      PlotLabel -> setup[[1, 1]]]; 
    \)

The flow plot tells you to expect all solutions to the sucked onto the t-axis.
Check it out:

\(Clear[y1, y2, y3, t]; 
    
{starter1, starter2, starter3} = {\(-3.0\), 2.0, 4.0}; 

    \(endtime = thigh; \)\)y1[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter1, 
            y[t], {t, 0, endtime}][[1
          ]]; y2[t_] = 
      y[t] /. NDSolve[diffeq /. starter -> starter2, 
            y[t], {t, 0, endtime}][[1
          ]]; \(y3[t_] = 
      y[t] /. \(NDSolve[diffeq /. starter -> starter3, 
            y[t], {t, 0, endtime}]\)[[1
          ]]; 


    solutionplots = 
      Plot[{y1[t], y2[t], y3[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        AspectRatio -> 1\/GoldenRatio, PlotRange -> {ylow, yhigh}, 
        AxesLabel -> {"t", "y[t]"}, 
        DisplayFunction -> Identity]; 


    Show[solutionplots, flowplot, 
      PlotLabel -> setup[[1, 1]], 
      DisplayFunction -> $DisplayFunction]; \)

If two solutions of this differential equation have nearby starting values on y[0], then the corresponding solutions will remain close to each other forever.
This differential equation is nearly insensitive to starter data.

T.2) The drinking versus driving model: Using a diffeq to analyze Bubba's toot

T.2.a)

Here are some facts:
-> A 12 - ounce long neck beer contains 13.6 grams of alcohol.
-> The typical human liver can eliminate 12 grams of alcohol per hour.
-> In all states, if your body contains more than 1 gram of alcohol per liter of body fluids, then you are too drunk to drive legally.
-> In some states, if your body contains more than 0.8 grams of alcohol per liter of body fluids, then you are too drunk to drive legally.                   
-> Bubba weighs 187 pounds.
-> A college-age male in good shape (as Bubba is) weighing K kilograms has about 0.68 K liters of fluids in his body.
To estimate the number of liters of fluids in Bubba's body, you use:

Needs["Miscellaneous`Units`"]Convert[187 PoundWeight, KilogramWeight]

Bubbafluids = 0.68 84.8217

This is in liters.

Here is Bubba's plan for the big party:  
He plans to drink plenty of beer at the party. But he is also worried about his driver's license. So he comes up with the scheme:
He'll drink beer gradually at the rate of one 12 - ounce long neck every 10 minutes for the first hour; then he'll nurse one 12 - ounce beer every half hour for the next two hours.  After all this beer, he'll stop drinking and just sit around looking cool until he's sure his blood alcohol concentration is down to a safe, legal level for driving.
Here is the function of time (in hours) that measures Bubba's intake rate of alcohol in grams per hour and a plot for the first 5 hours:

beer = 13.6; \(Clear[inrate, t]; 
inrate[t_] := 6\ beer /; 0 ≤ t < 1; 

    \(inrate[t_] := 2\ beer /; 1 ≤ t < 3; \)\)\(inrate[t_] := 0 /; 3 ≤ t; 


    Plot[inrate[t], {t, 0, 5}, 
      PlotStyle -> {{Thickness[0.01], YellowBrown}}, 
      PlotRange -> All, AspectRatio -> 1\/GoldenRatio, 
      AxesLabel -> {"t", "inrate[t]"}]; \)

Assuming Bubba's liver is normal, you can say the function of time (in hours) that measures his elimination rate of alcohol  in grams per hour is:

\(Clear[outrate, t]; 
outrate[t_] = 12; \)

Use the NDSolve instruction to help give a plot of the function that measures the alcohol (in grams) that will be found in Bubba's body anytime during the first 9 hours.
Then give a plot of the function that measures the concentration of alcohol in grams per liter in Bubba's body fluids.
Estimate how long Bubba will have to stay at the party before he becomes a legal driver.

Call the function that measures the alcohol (in grams) that will be found in Bubba's body anytime during the first 8 hours by the name 
       alcohol[t]. 
Assuming Bubba has no residue from the night before, you know that at the start of the party (t  =  0),
       alcohol[0]  =  0. 
You also know that 
       alcohol'[t]  =  inrate[t]  -  outrate[t].
Here's the plot of alcohol[t] for the first 9 hours:

beer = 13.6; \(Clear[inrate, outrate, t]; 

    inrate[t_] := 6\ beer /; 0 ≤ t < 1; 

    \(inrate[t_] := 2\ beer /; 1 ≤ t < 3; \)\)inrate[t_] := 0 /; 3 ≤ t; outrate[t_] = 12; a = 0; b = 9; starter = 0;

Clear[solution, t, y, alcohol];
solution=NDSolve[{y'[t]==inrate[t] - outrate[t],y[0] == starter},y[t],{t, a, b}];
alcohol[t_] = 
        y[t] /. solution[[1]];

Plot[alcohol[t], {t, a, b}, 
        PlotStyle -> {{YellowBrown, Thickness[0.01]}}, 
        AxesLabel -> {"t", "alcohol[t]"}, 
        AxesOrigin -> {a, starter}, AspectRatio -> 1, 
        PlotRange -> All];

Call the function that measures the concentration of alcohol in grams per liter in Bubba's body fluids by the name concentration[t]. This function is given by:

\(Clear[concentration]; 
concentration[t_] = alcohol[t]\/Bubbafluids\)

Here comes the plot showing concentration[t] together with the legal thresholds 1.0 grams per liter, 0.8 grams per liter, and the conservative threshold 0.5 grams per liter:

\(\(Plot[{concentration[t], 1.0, 0.8, 0.5}, {t, a, b}, 
      PlotStyle -> {{YellowBrown, Thickness[0.01]}, Red, Red, Green}, 

      	AspectRatio -> 1, 
      AxesLabel -> {"t", "concentration[t]"}, 
      PlotRange -> All]; \)\)

In some states, Bubba can legally drive home about 7 hours after the party begins.  In other states, he must stay at the party for nearly eight hours. 
Just to be on the safe side, it might be a good idea for Bubba to stay at the party a full nine hours.
Go back and run a drinking schedule you or a friend is likely to use.  You might pick up some useful information.

T.3) Using parametric plotting to plot the predator population as a function of the prey population

Remember the predator-prey model from the previous lesson.  In case you've forgotten about it, here's the scoop again:

Two species coexist in a closed environment. One species, the predator, feeds on the other, the prey. There is always plenty of food for the prey, but the predators eat nothing but the prey.
Put 
            pred[t]  =  population  of  predators  at  time  t.
Put  
            prey[t]  =  population  of  prey  at  time  t.

It's reasonable to assume that there are positive constants a and b such that:
           prey'[t]  =  a prey[t]  -  b prey[t] pred[t]
because the abundance of food for the prey allows the birth rate of the prey to be proportional to their current number, and the death rate of prey is proportional to both the current number of prey and the current number of predators.

It also makes some sense to assume that there are positive constants c and d such that:
             pred'[t]  =  -c pred[t]  +  d pred[t] prey[t]
because it's reasonable to assume that the death rate of the predators is likely to be proportional to the current population of predators and that the birth rate of the predators is proportional to both the current number of the predators and the size of the food supply (the prey).

Here is what happens for a sample choice of a,  b,  c,  d with the prey starting off with a population of 4 units and the predators starting out with a population of 1 unit.

Here is Mathematica's fake plot of the prey population as a function of time[t] for the first 50 time units:

endtime = 50; a = 0.7; b = 0.3; c = 0.3; d = 0.1;

Clear[pred, prey, t];
approxsolutions=NDSolve[{prey'[t]==a prey[t] - b prey[t] pred[t],pred'[t]==-c pred[t] + d pred[t] prey[t],prey[0] == 4,pred[0] == 1},{prey[t], pred[t]},{t, 0, endtime}];

Clear[fakepred, fakeprey];
fakepred[t_] = 
        pred[t] /. 
          approxsolutions[[1]];
fakeprey[t_] = 
        prey[t] /. 
          approxsolutions[[1]];

preyplot = 
        Plot[fakeprey[t], {t, 0, endtime}, 
          PlotStyle -> {{Blue, Thickness[0.015]}}, 
          AxesLabel -> {"t", "prey"}, AspectRatio -> (1/3)]
          ;

And the predator population as a function of time t:

predatorplot = 
      Plot[fakepred[t], {t, 0, endtime}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"t", "predators"}, 
        AspectRatio -> (1/3)];

Here they are together:

both = Show[predatorplot, preyplot];

Note the relationship between the predator population's crests and dips and prey population's crests and dips.  Just after the predators start growing, the prey is just about eaten up.

T.3.a)

Use parametric plotting to display the predator population as a function of the prey population and describe what you see.

Here it is:

fullplot = 
      ParametricPlot[{fakeprey[t], fakepred[t]}, {t, 0, endtime}, 
        PlotStyle -> {{Thickness[0.01], DarkGreen}}, 
        PlotRange -> All, 
        AxesLabel -> {"prey", "predator"}];

A distorted ellipse that some folks call a closed curve.
This is a little hard to interpret. 
To get a grasp on it, see what happens as t goes from 0 to 3:

three = 
      ParametricPlot[{fakeprey[t], fakepred[t]}, {t, 0, 3}, 
        PlotStyle -> {{Thickness[0.01], DarkGreen}}, 
        AxesLabel -> {"prey", "predator"}, 
        PlotRange -> {{0, 9}, {0, 5}}, PlotLabel -> 3];

Now see what happens as t goes from 0 to 7:

seven = 
      ParametricPlot[{fakeprey[t], fakepred[t]}, {t, 0, 7}, 
        PlotStyle -> {{Thickness[0.01], DarkGreen}}, 
        AxesLabel -> {"prey", "predator"}, 
        PlotRange -> {{0, 9}, {0, 5}}, PlotLabel -> 7];

Now see what happens as t goes from 0 to 12:

twelve = 
      ParametricPlot[{fakeprey[t], fakepred[t]}, {t, 0, 12}, 
        PlotStyle -> {{Thickness[0.01], DarkGreen}}, 
        AxesLabel -> {"prey", "predator"}, 
        PlotRange -> {{0, 9}, {0, 5}}, PlotLabel -> 12];

Now see what happens as t goes from 0 to 15:

fifteen = 
      ParametricPlot[{fakeprey[t], fakepred[t]}, {t, 0, 15}, 
        PlotStyle -> {{Thickness[0.01], DarkGreen}}, 
        AxesLabel -> {"prey", "predator"}, 
        PlotRange -> {{0, 9}, {0, 5}}, PlotLabel -> 15];

Holy smoke!
This is the same as:

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}];

The above plot shows what happens as t goes from 0 to 50.
Review the plots:

Show[GraphicsArray[{{three, seven}, {twelve, fifteen}}]];

As time advances, the parametric points {prey[t], pred[t]} advance around the closed curve in a counterclockwise way.  As time advances on and on, the parametric points {prey[t], pred[t]} cycle around this curve over and over again. 
Take another look but this time setting the axes' origin to {(a/b), (c/d)}:

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}, 
      AxesOrigin -> {(c/d), (a/b)}];

Again the closed curve tells you that predator and prey populations are periodic; they go through repeated cycles. And the time cycle for each is the same. The four sectors above defined by the axes indicate four phases depicting trends that reverse themselves as the curve crosses the axes at prey = (c/d) and predator = (a/b).

T.3.b.i)

Wait a minute!
Take another look at the last plot.

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}, 
      AxesOrigin -> {(c/d), (a/b)}];

Those axes pierce the curve at the maximum and minimum values of prey and predator populations.
Is this just an accident?

Stupid question.

T.3.b.ii)

Explain why those axes pierce the plot at the maximum and minimum values of prey and predator populations.

Look at the original differential equations
        prey'[t]  =    a prey[t]  -  b prey[t] pred[t],
        pred'[t]  =  -c pred[t]  +  d pred[t] prey[t].
Take the first:
        prey'[t]  =    a prey[t]  -  b prey[t] pred[t].
At the times t at which the prey population is at its maximum or at its minimum, you gotta have 
        prey'[t]  =  0.
Consequently, at the times t at which the prey population is at its maximum or at its minimum,
          0  =  a prey[t]  -  b prey[t] pred[t],
              =  prey[t] (a   -  b pred[t]).
This tells you that at the times t at which the prey population is at its maximum or at its minimum,
             pred[t] = (a/b).

Similarly, use
           pred'[t]  =  -c pred[t]  +  d pred[t] prey[t] 
to see that for a maximum or minimum predator population, you have
           0  =  -c pred[t]  +  d pred[t] prey[t]
               =  pred[t] (-c   +  d  prey[t]),
so when the predator population is largest or smallest,
           prey[t] = (c/d).

This explains why when you use AxesOrigin -> {(c/d), (a/b)}; the  axes pierce the plot at the maximum and minimum values of prey and predator populations.

Show[fullplot, PlotRange -> {{0, 9}, {0, 5}}, 
      AxesOrigin -> {(c/d), (a/b)}];

Math happens again.

T.4) More on reading the predator-prey model

T.4.a)

Here is what happens in the predator-prey model for the sample choice of proportionality constants 
         a  =  0.2,  b  =  0.05,  c  =  0.28,  and  d = 0.04
with the prey starting off with a population of 20 units and the predators starting out with a population of 2 units.

\(endtime\  = \ 50; 

a\  = \ 0.2; 
b\  = \ 0.05; 
c\  = \ 0.28; 

    d\  = \ 0.04; 

Clear[pred, prey, t]; 

    approxsolutions\  = \ 

      NDSolve[{\(prey'\)[t]\  == \ \ a\ prey[t]\  - \ b\ prey[t]\ pred[t], 

          \ \ \ \ \ \ \ \ \ 
          \(pred'\)[t]\  == \ \(-c\)\ pred[t]\  + \ d\ pred[t]\ prey[t], 

          \ \ \ \ \ \ \ \ \ prey[0]\  == \ 20, 
\ \ \ \ \ \ \ \ \ 
          pred[0]\  == \ 2}, 
\ \ \ {prey[t], pred[t]}, {t, 0, endtime}]; 

    \ \ \ 
Clear[fakepred, fakeprey]; 

    fakepred[t_]\  = \ pred[t] /. approxsolutions[\([1]\)]; 


    fakeprey[t_]\  = \ prey[t] /. approxsolutions[\([1]\)]; 


    predpreyplot\  = 

      Plot[{fakeprey[t], fakepred[t]}, {t, 0, endtime}, 
\ \ \ \ \ \ \ 
        PlotStyle -> {{Blue, Thickness[0.015]}, 

            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Red, Thickness[0.007]}}, 

        \ \ \ \ \ \ \ AxesLabel -> {"t", ""}]; \)

Add the horizontal line that hits the vertical axis at (a/b) to this plot. Inspect the plot; then describe what you see and explain why you see it.

Here you go:

Show[predpreyplot, Graphics[Line[{{0, (a/b)}, {endtime, (a/b)}}]]];

Interesting. 
When pred[t]  > (a/b), then prey[t] goes down.
When pred[t]  < (a/b), then prey[t] goes up. 
And prey[t] takes its largest and smallest values when
         pred[t]  =  (a/b).

To explain this phenomenom, look at one of the differential equations that went into this plot:
       prey'[t]  =   a prey[t]  -  b prey[t] pred[t]
                         =   prey[t] (a  -  b  pred[t]).
When pred[t]  > (a/b), then (a  -  b  pred[t])  <  0.
So when pred[t]  > (a/b), then 
       prey'[t]  =   prey[t] (a  -  b  pred[t])  <  0. 
That's why the prey curve goes down when pred[t]  >  (a/b).

Similarly when pred[t]  <  (a/b), then prey'[t]  >  0. 
That's why the prey curve goes up when pred[t]  < (a/b).
When pred[t]  = (a/b), then
          prey'[t]  =   
        prey[t] (a  -  b  pred[t])  =  0.
That's why prey[t] takes its largest and smallest values when
       pred[t]  =  (a/b).

Amazing what a little calculus can do.