The phase plane

Eigenvalues, eigenvectors, and equilibrium points

- Marcus McGuff, Austin Community College

Manipulate[

Refresh[If[var===Null,Clear[var];var["highrect"]={};var["computePlotrange"]=plotrange;
var["xEqOld"]=xEq;
var["yEqOld"]=yEq;
var["paramsOld"]=params];
(*This is the list of functions to be included in the drop-down lists; can add extra info later inside each inner bracket; be SURE that each function is enclosed in ToBoxes[...], and function is entered in StandardForm
*)
If[xText===Null,Clear[x,y];
xFunctions={
{ToBoxes[a x + b y]},
{ToBoxes[x(a - b y)]},
{ToBoxes[x(1 - a x - b y)]},
{ToBoxes[a y]}
};
xText=xFunctions[[1,1]]];
If[yText===Null,Clear[x,y];
yFunctions={
{ToBoxes[c x + d y]},
{ToBoxes[y(c x - d)]},
{ToBoxes[y(1 - c y - d x)]},
{ToBoxes[c x - d (x^3)]}
};
yText=yFunctions[[1,1]]],None];

(* Convert input text boxes into expressions in TraditionalForm*)
xEq=ToExpression[xText,TraditionalForm];
yEq=ToExpression[yText,TraditionalForm];
If[xEq===Null&&yEq===Null,var["xEqOld"]=xEq;
var["yEqOld"]=yEq];


(* Finds parameters in xEq and yEq, then writes the replacement rules for them; can specify exclusion cases here; set default value for all parameters here with {_->1.} *)
params=SortBy[Thread[Rule[#,Replace[#,Join[params,{_->1.}],1]]],First]&[Cases[Level[{xEq,yEq},-1],s_Symbol/;And[FreeQ[{x,y},s],FreeQ[Attributes[s],Locked|Constant]]]];
If[params===Null,var["paramsOld"]=params];

(* Compute the critical points and Eignesets *)
Refresh[
If[findCrPtsYN&&PolynomialQ[xEq,{x,y}]&&PolynomialQ[yEq,{x,y}],
eqPts={};
var["allEqPts"]=Union[Quiet@Solve[Thread[{xEq,yEq}=={0,0}]/.params,{x,y}]];

If[!NumberQ[Simplify[(xEq/yEq)/.params]],
Table[If[({x,y}/.var["allEqPts"][[k]])∈Reals,AppendTo[eqPts,var["allEqPts"][[k]]]],{k,1,Length[var["allEqPts"]]}],
AppendTo[eqPts,{x->0,y->0}]];

If[Length[eqPts]>0,
localLin=D[{xEq,yEq}/.params,{{x,y}}]/.eqPts//N;
eigenSets=Map[Eigensystem,localLin],
findCrPtsYN=False;MessageDialog["There are no Real-valued critical points."]],
If[findCrPtsYN,findCrPtsYN=False;MessageDialog["Your equations need to be polynomials in x and y to find the critical points in this demonstration."]]],TrackedSymbols:>{xText,yText,params}];


If[!(var["xEqOld"]===xEq&&var["yEqOld"]===yEq&&var["paramsOld"]===params),

(* Recompute existing trajectories whenever params change *)
If[Length[pts]>0,
Table[
pts[[k,2]]=First[({x[t],y[t]}/.Quiet[var["eqSol"]=NDSolve[{Thread[{x'[t],y'[t]}==({xEq,yEq}/.params/.{x->x[t],y->y[t]})],Thread[{x[0],y[0]}==pts[[k,1]]]},{x,y},{t,0,tMax}]])];
pts[[k,3]]=Flatten[Thread@DifferentialEquations`InterpolatingFunctionAnatomy`InterpolatingFunctionDomain[Flatten[{x,y}/.var["eqSol"]]],1],
{k,1,Length[pts]}]];

var["xEqOld"]=xEq;
var["yEqOld"]=yEq;
var["paramsOld"]=params];



EventHandler[Grid[{{Text@Style[If[findCrPtsYN&&(Length[eqPts]>0),Row[{"Equilibrium pts: ",Dynamic[{x,y}/.eqPts]}],"Phase plane"]]},
{Tooltip[Dynamic[Graphics[{
(* Vector field or StreamPlot *)
Dynamic[First[plotter[Evaluate[{xEq,yEq}/.params],
{x,Part[var["computePlotrange"],1,1],Part[var["computePlotrange"],1,2]},{y,Part[var["computePlotrange"],2,1],Part[var["computePlotrange"],2,2]},
StreamStyle->Red,
VectorScale->{0.03,Automatic,1&},
VectorStyle->{Blue,Arrowheads[0.02]}]],
(*,*)
(*TrackedSymbols->{pts,plotter,xText,yText,var["computePlotrange"]}*)],

(* Set of solution curves and initial points, with Tooltips displaying x(t) and y(t) *)
Dynamic[If[Length[pts]>0,
Table[
Tooltip[{{Red,PointSize[Medium],Point[pts[[k,1]]]},
First@ParametricPlot[pts[[k,2]],{t,Max[pts[[k,3,All,1]]],Min[pts[[k,3,All,2]]]}]},
With[{k=k},Dynamic@Column[
{Plot[Part[pts,k,2,1],{t,pts[[k,3,1,1]],pts[[k,3,1,2]]},AxesLabel->{t,x},AxesOrigin->{0,0},
PlotRange->2{Part[var["computePlotrange"],1,1],Part[var["computePlotrange"],1,2]}],
Plot[Part[pts,k,2,2],{t,pts[[k,3,2,1]],pts[[k,3,2,2]]},AxesLabel->{t,y},AxesOrigin->{0,0},
PlotRange->2{Part[var["computePlotrange"],2,1],Part[var["computePlotrange"],2,2]}]}]]],
{k,1,Length[pts]}],
{}]],

(* Highlight rectangle for zooming *)
Dynamic[var["highrect"]],

(* Plot equilibrium points and draw eigenvectors *)
Dynamic[If[findCrPtsYN&&(Length[eqPts]>0),
{PointSize[Large],Red,Thickness[0.004],Dashed,
Table[With[{ptk={x,y}/.eqPts[[k]]},
Append[If[(eigenSets[[k,1,1]]∈Reals),
Table[Tooltip[Line[{ptk-10eigenSets[[k,2,m]],ptk+10eigenSets[[k,2,m]]}],
Text@Style[Row[{"λ = ",eigenSets[[k,1,m]],", Eigenvector: ",eigenSets[[k,2,m]]}],Medium]],{m,1,2}],{}],
Tooltip[Point[ptk],With[{k=k},Dynamic@Column[{Text@Style[Row[{"λ = ",eigenSets[[k,1]]}],Medium],StreamPlot[localLin[[k]].{x,y},{x,-0.25(Part[var["computePlotrange"],1,2]-Part[var["computePlotrange"],1,1]),0.25(Part[var["computePlotrange"],1,2]-Part[var["computePlotrange"],1,1])},{y,-0.25(Part[var["computePlotrange"],2,2]-Part[var["computePlotrange"],2,1]),0.25(Part[var["computePlotrange"],1,2]-Part[var["computePlotrange"],1,1])},
Axes->False,Frame->False,PerformanceGoal->"Speed"]}]]]]],
{k,1,Length[eqPts]}]},
{}]]

},
PlotRange->Dynamic[plotrange],AspectRatio->1,AxesOrigin->{0,0},Axes->True,AxesLabel->{x,y},Frame->False,ImageSize-> 500,FilterRules[Options[VectorPlot],Options[Graphics]]
]],
Dynamic[Text@Style[MousePosition["Graphics"],Medium]]]},
{Text@Grid[{{"Mouse actions",SpanFromLeft},{"Click - Plot solution curve","","Shift/Drag - Move graph"},
{"Alt(Cmd)/Drag - Zoom in","","Shift/Alt(Cmd)/Drag - Zoom out"}},Alignment->Center,Dividers->{False,{False,True,False,False}},ItemStyle->Directive[Black,10,FontFamily->"Times"]]}}],

{"MouseClicked":>
Switch[mouseMode[],
"pan",None,
"zoomIn",None,
"zoomOut",None,
"special",If[!MemberQ[pts[[All,1]],MousePosition["Graphics"]],AppendTo[pts,
{MousePosition["Graphics"],First[({x[t],y[t]}/.Quiet[var["eqSol"]=NDSolve[{Thread[{x'[t],y'[t]}==({xEq,yEq}/.params/.{x->x[t],y->y[t]})],Thread[{x[0],y[0]}==MousePosition["Graphics"]]},{x,y},{t,0,tMax}]])],
Flatten[Thread@InterpolatingFunctionDomain[Flatten[{x,y}/.var["eqSol"]]],1]}],
None]],
"MouseDown":>
(readMouse["start"];
Switch[mouseMode[],
"pan",None,
"zoomIn",None,
"zoomOut",None,
"special",None]
),
"MouseDragged":>
(readMouse["end"];
Switch[mouseMode[],
"pan",(plotrange=plotrange-Apply[Subtract,-plotrange,1](var["endScaled"]-var["startScaled"]);
readMouse["start"]),
"zoomIn",var["highrect"]={Opacity[0.2],Blue,Rectangle[var["start"],var["end"]]},
"zoomOut",var["highrect"]={Opacity[0.2],Red,Rectangle[var["start"],var["end"]]},
"special",None]),
"MouseUp":>(Switch[mouseMode[],
"pan",None,
"zoomIn",If[var["highrect"]=!={},
plotrange=Sort/@Transpose[var/@{"start","end"}]],
"zoomOut",If[var["highrect"]=!={},
plotrange=1/Abs[(var["endScaled"]-var["startScaled"])]plotrange],
"special",None];
var["highrect"]={};
readMouse["start"];
pts=Union@pts;
var["computePlotrange"]=plotrange)}],

(* Allows for TraditionalForm input - Must apply ToExpression[xText,TraditionalForm] before use, see xEq and yEq below *)
OpenerView[{
Row[{Control[{{xText,,"(dx/dt)="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@xFunctions[[All,1]],ControlType->PopupMenu}],"                                    ",
Control[{{yText,,"(dy/dt)="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@yFunctions[[All,1]],ControlType->PopupMenu}]}],
If[$LicenseID≠"L9999-9999",
Row[{Control[{{xText,,"(dx/dt)="},TraditionalForm[InputField[Dynamic[xText],Boxes]]&}],"  ",Control[{{yText,,"(dy/dt)="},TraditionalForm[InputField[Dynamic[yText],Boxes]]&}]}],Text@Style[
"Extra input options not available in the free Mathematica Player; to enter your own equations, you must be using a full version of Mathematica",Italic]]}],
(*OpenerView[{
Row[{Control[{{xEq,,"(dx/dt)="},xFunctions[[All,1]],ControlType->PopupMenu}],"                             ",Control[{{yEq,,"(dy/dt)="},yFunctions[[All,1]],ControlType->PopupMenu}]}],
Row[{Control[{{xEq,,"(dx/dt)="},ControlType->InputField}],"    ",Control[{{yEq,,"(dy/dt)="},ControlType->InputField}]}]}],*)

(* Allows tMax to continue moving upwards indefinitely; tLimit is initial upper value for tMax *)
{{tLimit, 20}, ControlType -> None}, {{tMax, 5,"Max time domain"}, Refresh[Slider[Dynamic[tMax, {Automatic, (tLimit = tMax + 0.5tLimit) &}], {0.000000001, tLimit},Appearance->"Labeled"], TrackedSymbols -> {tLimit}] &}, 

Row[{Button["Clear Curves",pts={},BaseStyle->Small,ImageSize->All],"     ",
Button["Default Range",plotrange={{-5,5},{-5,5}},BaseStyle->Small,ImageSize->All],"     ",
Control@{{findCrPtsYN,False,"Show equilibrium pts"},{True,False}},"   ",Control@{{plotter,VectorPlot,"Show me all"},{VectorPlot,StreamPlot},Checkbox}
}],

Delimiter,
(* Produces sliders dynamically for parameters *)
Item[Column[{"Parameters",Dynamic[Grid[MapIndexed[{params[[First@#2,1]],Manipulator[Dynamic[params[[First@#2,2]]],{-5,5,0.01},Appearance->{"Labeled","Open"}]}&,params]]]}],ControlPlacement->Left],

{xEq,ControlType->None},
{yEq,ControlType->None},
{{params,{}},ControlType->None},
{{plotrange,{{-5,5},{-5,5}}},ControlType->None},

(* temporary holding variable *)
{var,ControlType->None},

(* pts -> list of info for each init cond: { coord , {x[t],y[t]}, {{t domain for x}, {t domain for y}}
so: pts[[k,1]] gives initial condition, pts[[k,2]] gives solution curves, and pts[[k,3]] gives safe domains of t for each solution curve*)
{{pts,{}},ControlType->None},

{{eqPts,{}},ControlType->None},
{{localLin,{}},ControlType->None},
{{eigenSets,{}},ControlType->None},
{xFunctions,None},
{yFunctions,None},

Initialization:>(mouseMode[]:=Switch[CurrentValue/@{"ShiftKey","CommandKey"},
{False,False},"special",
{True,False},"pan",
{False,True},"zoomIn",
{True,True},"zoomOut"];
readMouse[startOrEnd_]:=(var[startOrEnd]=MousePosition["Graphics"];
var[startOrEnd<>"Scaled"]=MousePosition["GraphicsScaled"]);
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]
),Alignment->Center,SynchronousInitialization->False]