Systems in 3 dimensions

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,z];
xFunctions={
{ToBoxes[x+y-z]},
{ToBoxes[-x+y-z]},
{ToBoxes[a x + b y+c z]},
{ToBoxes[s (-x+y)]}
};
xText=xFunctions[[1,1]]];
If[yText===Null,Clear[x,y,z];
yFunctions={
{ToBoxes[2y+3z]},
{ToBoxes[-2y+3z]},
{ToBoxes[p x + q y+r z]},
{ToBoxes[r x-y-x z]}
};
yText=yFunctions[[1,1]]],None];
If[zText===Null,Clear[x,y,z];
zFunctions={
{ToBoxes[3z]},
{ToBoxes[-3z]},
{ToBoxes[3y-3z]},
{ToBoxes[-2y+2z]},
{ToBoxes[u x + v y+w z]},
{ToBoxes[-b z+x y]}
};
zText=zFunctions[[1,1]]];

plotrange={{-(zoom^3),(zoom^3)},{-(zoom^3),(zoom^3)},{-(zoom^3),(zoom^3)}};

(* Convert input text boxes into expressions in TraditionalForm*)
xEq=ToExpression[xText,TraditionalForm];
yEq=ToExpression[yText,TraditionalForm];
zEq=ToExpression[zText,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,zEq},-1],s_Symbol/;And[FreeQ[{x,y,z},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,z}]&&PolynomialQ[yEq,{x,y,z}]&&PolynomialQ[zEq,{x,y,z}],
eqPts={};
var["allEqPts"]=Union[Quiet@Solve[Thread[{xEq,yEq,zEq}=={0,0,0}]/.params,{x,y,z}]];

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

If[Length[eqPts]>0,
localLin=D[{xEq,yEq,zEq}/.params,{{x,y,z}}]/.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,zText,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],z[t]}/.Quiet[var["eqSol"]=NDSolve[{Thread[{x'[t],y'[t],z'[t]}==({xEq,yEq,zEq}/.params/.{x->x[t],y->y[t],z->z[t]})],Thread[{x[0],y[0],z[0]}==pts[[k,1]]]},{x,y,z},{t,0,tMax}]])];
pts[[k,3]]=Flatten[Thread@DifferentialEquations`InterpolatingFunctionAnatomy`InterpolatingFunctionDomain[Flatten[{x,y,z}/.var["eqSol"]]],1],
{k,1,Length[pts]}]];

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

curPt={Append[xyPt,zPt],First[({x[t],y[t],z[t]}/.Quiet[var["eqSol"]=NDSolve[{Thread[{x'[t],y'[t],z'[t]}==({xEq,yEq,zEq}/.params/.{x->x[t],y->y[t],z->z[t]})],Thread[{x[0],y[0],z[0]}==Append[xyPt,zPt]]},{x,y,z},{t,0,tMax}]])],
Flatten[Thread@DifferentialEquations`InterpolatingFunctionAnatomy`InterpolatingFunctionDomain[Flatten[{x,y,z}/.var["eqSol"]]],1]};


Grid[{{Text@Style[Grid[{If[findCrPtsYN&&(Length[eqPts]>0),List["Equilibrium pts: ",Dynamic[{x,y,z}/.eqPts]],List["Phase space",SpanFromLeft]],{"Initial condition: ",Dynamic@Append[xyPt,zPt]}}]]},
{Dynamic[Graphics3D[{
(* Vector field or StreamPlot *)
(*Dynamic[First[VectorPlot3D[Evaluate[{xEq,yEq,zEq}/.params],
{x,Part[plotrange,1,1],Part[plotrange,1,2]},{y,Part[plotrange,2,1],Part[plotrange,2,2]},
{z,Part[plotrange,3,1],Part[plotrange,3,2]},
VectorPoints->5,VectorScale->Medium,
VectorStyle->Blue]]
],
*)

Tooltip[{{Red,PointSize[Large],Point[curPt[[1]]]},
First@ParametricPlot3D[curPt[[2]],{t,Max[curPt[[3,All,1]]],Min[curPt[[3,All,2]]]},PlotStyle->Thick]},
Dynamic@Column[
{Plot[Part[curPt,2,1],{t,curPt[[3,1,1]],curPt[[3,1,2]]},AxesLabel->{t,x},AxesOrigin->{0,0},
PlotRange->2{Part[plotrange,1,1],Part[plotrange,1,2]}],
Plot[Part[curPt,2,2],{t,curPt[[3,2,1]],curPt[[3,2,2]]},AxesLabel->{t,y},AxesOrigin->{0,0},
PlotRange->2{Part[plotrange,2,1],Part[plotrange,2,2]}],
Plot[Part[curPt,2,3],{t,curPt[[3,3,1]],curPt[[3,3,2]]},AxesLabel->{t,z},AxesOrigin->{0,0},
PlotRange->2{Part[plotrange,3,1],Part[plotrange,3,2]}]}]],


(* Set of solution curves and initial points, with Tooltips displaying x(t) and y(t) and z(t) *)
Dynamic[If[Length[pts]>0,
Table[
Tooltip[{{Red,PointSize[Medium],Point[pts[[k,1]]]},
First@ParametricPlot3D[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[plotrange,1,1],Part[plotrange,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[plotrange,2,1],Part[plotrange,2,2]}],
Plot[Part[pts,k,2,3],{t,pts[[k,3,3,1]],pts[[k,3,3,2]]},AxesLabel->{t,z},AxesOrigin->{0,0},
PlotRange->2{Part[plotrange,3,1],Part[plotrange,3,2]}]}]]],
{k,1,Length[pts]}],
{}]],


(* Show nullclines *)
Dynamic[If[showNullclinesYN,
First@ContourPlot3D[Evaluate[Thread[{xEq,yEq,zEq}=={0,0,0}]/.params],{x,Part[plotrange,1,1],Part[plotrange,1,2]},{y,Part[plotrange,2,1],Part[plotrange,2,2]},
{z,Part[plotrange,3,1],Part[plotrange,3,2]},MaxRecursion->ControlActive[3,Automatic],Mesh->None,
ContourStyle->
{Directive[Green,Opacity[0.3]],
Directive[Orange,Opacity[0.3]],
Directive[Yellow,Opacity[0.3]]}],
{}]],

(* Plot equilibrium points and draw eigenvectors *)
Dynamic[If[findCrPtsYN&&(Length[eqPts]>0),
{PointSize[Large],Blue,Thickness[0.003],Dashed,
Table[With[{ptk={x,y,z}/.eqPts[[k]]},
Append[
Append[Table[If[(eigenSets[[k,1,m]]∈Reals),
Tooltip[Line[{ptk-(zoom^3)eigenSets[[k,2,m]],ptk+(zoom^3)eigenSets[[k,2,m]]}],
Text@Style[Row[{"λ = ",eigenSets[[k,1,m]],", Eigenvector: ",eigenSets[[k,2,m]]}],Medium]],{}],{m,1,Length[eigenSets[[k,1]]]}],Red],

Tooltip[Point[ptk],With[{k=k},Dynamic@Column[{Text@Style[Row[{"λ = ",eigenSets[[k,1]]}],Medium]}]]]]],
{k,1,Length[eqPts]}]},
{}]]

},
PlotRange->Dynamic[plotrange],AxesOrigin->{0,0,0},AxesLabel->{x,y,z},ImageSize-> 500,FilterRules[Options[VectorPlot3D],Options[Graphics3D]]
]]},

{Text@Grid[{
If[findCrPtsYN,{"Red - Eq. pts/initial conditions","","Blue - Eigenvectors"},{}],
If[showNullclinesYN,{"Green - x nullclines","","Orange - y nullclines","","Yellow - z nullclines"},{}]},Alignment->Center,Dividers->{False,{True,False,False,False}},ItemStyle->Directive[Black,10,FontFamily->"Times",Italic]]}}],


(* 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}],"  ",
Control[{{zText,,"(dz/dt)="},(#->ToExpression[#,StandardForm,TraditionalForm])&/@zFunctions[[All,1]],ControlType->PopupMenu}]}],
Dynamic@If[StringCases[$InstallationDirectory,"Player"]=={},
Row[{Control[{{xText,,"(dx/dt)="},TraditionalForm[InputField[Dynamic[xText],Boxes]]&}],"  ",Control[{{yText,,"(dy/dt)="},TraditionalForm[InputField[Dynamic[yText],Boxes]]&}],"  ",Control[{{zText,,"(dz/dt)="},TraditionalForm[InputField[Dynamic[zText],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]]}],

(* 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}] &}, 


(*Grid[{{Text@Style["# of curves in each direction"],SpanFromLeft,SpanFromLeft},
{Control@{{xGrid,1,"X"},1,10,1},
Control@{{yGrid,1,"Y"},1,10,1},
Control@{{zGrid,1,"Z"},1,10,1}}}],

{x0,ControlType->None},
{y0,ControlType->None},
{z0,ControlType->None},
*)

Row[{"     ",
Control@{{findCrPtsYN,False,"Show equilibrium pts"},{True,False}},"   ",
Control@{{showNullclinesYN,False,"Show nullclines"},{True,False}}
}],

Item[Column[{Button["Clear Curves",pts={},BaseStyle->Small,ImageSize->All],"     ",
Button["Add curve",
If[!MemberQ[pts[[All,1]],Append[xyPt,zPt]],AppendTo[pts,
{Append[xyPt,zPt],First[({x[t],y[t],z[t]}/.Quiet[var["eqSol"]=NDSolve[{Thread[{x'[t],y'[t],z'[t]}==({xEq,yEq,zEq}/.params/.{x->x[t],y->y[t],z->z[t]})],Thread[{x[0],y[0],z[0]}==Append[xyPt,zPt]]},{x,y,z},{t,0,tMax}]])],
Flatten[Thread@InterpolatingFunctionDomain[Flatten[{x,y,z}/.var["eqSol"]]],1]}],
None],
BaseStyle->Small,ImageSize->All],
Text@Style["Initial conditions:"],Text@Style["Horizontal"],
Dynamic@Control@{{xyPt,{1,1},""},{-(zoom^3),-(zoom^3)},{(zoom^3),(zoom^3)},ControlPlacement->Right},
Grid[{{Text@Style["Vertical"],Text@Style["Zoom"]},
{Dynamic@Control@{{zPt,1,""},-(zoom^3),(zoom^3),ControlType->VerticalSlider},Control@{{zoom,2,""},0.1,10,ControlType->VerticalSlider}}}]}],ControlPlacement->Right],

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},
{zEq,ControlType->None},
{{params,{}},ControlType->None},
{{plotrange,{}},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},
{curPt,ControlType->None},

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

Initialization:>(computeIC[start_,end_,totalDivs_,index_]:=start+(index-0.5)(end-start)/totalDivs;
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]