NullclinePlot
In[]:=
Options[NullclinePlot] = { NullclineContourPlotPoints->50, NullclineShadingMethod -> Automatic, ContourStyle->Automatic, NullclineShading->False, ShowEquilibria->False, EquilibriumPointColor->GrayLevel[0], Nullclines->True, NullclineStyle->{{},{}}}; TurnPts[set_] := First /@ set[[(1+(Flatten[Position[Apply[Times, Partition[DeleteCases[Sign[Apply[#1 - #2 &, Partition[First /@ set, 2, 1], 2]], 0], 2, 1], 2], -1]]))]] /; Depth[set] == 3 TurnPts[sets_] := Apply[Join, TurnPts /@ sets] /; Depth[sets] == 4 ToFunctions[set_] := If[#[[1,1]] > #[[-1,1]], Reverse[#], #]& /@ (set[[#]] & /@ Apply[Range, Partition[Prepend[Append[1+(Flatten[Position[Apply[Times,Partition[ DeleteCases[ Sign[Apply[#1 - #2 &, Partition[First /@ set, 2, 1], 2]], 0], 2,1], 2], -1]]), Length[set]], 1], 2, 1], 2]); BreakAtXVal[set_, x_?NumberQ] := Module[{n = Length[set], i = home[First[Transpose[set]], x]}, If[1 <= i < n, linInt = {{x, set[[i, 2]] + (x - set[[i,1]])* (set[[i+1, 2]]-set[[i,2]])/(set[[i+1, 1]] - set[[i,1]])}}; {Join[Take[set, i], If[i<=n-1, linInt, {}]], Join[If[i>0, linInt, {}], Take[set, i-n]]}, {set}] ] /; Depth[set] < 4 BreakAtXVal[sets_, x_?NumberQ] := (Apply[Join, BreakAtXVal[#, x] & /@ sets]) /; Depth[sets] == 4 BreakAtXVal[set_, xVals_List] := Fold[ DeleteCases[BreakAtXVal[#1, #2], {}]&, set, xVals] /; Depth[set] == 3 BreakAtXVal[sets_, xVals_List] := (Apply[Join, BreakAtXVal[#, xVals] & /@ sets]) /; Depth[sets] == 4 home[xvals_, l_] := Infinity /; l >= Last[xvals]; home[xvals_, l_] := -Infinity /; l <= First[xvals]; home[{r_, s_}, _] := 1 /; r < s; home[xvals_, x_] := Module[{n = Length[xvals], h}, h = Quotient[n, 2]; Which[x < xvals[[h]], home[Take[xvals, h], x], x > xvals[[h+1]], h + home[Take[xvals, h - n], x], True, h] ] home1[{r_, s_}, _] := 1 /; r < s; home1[xvals_, x_] := Module[{n = Length[xvals], h}, h = Quotient[n, 2]; Which[x < xvals[[h]], home1[Take[xvals, h], x], x > xvals[[h+1]], h + home1[Take[xvals, h - n], x], True, h] ] NullclinePlot[funcs_List,{x_,xmin_,xmax_},{y_,ymin_,ymax_},opts___]:= Module[{}, {cpp,shadedQ,stabQ, cs, stabCol, showLines, shadmethod, nullcSty} = {NullclineContourPlotPoints, NullclineShading, ShowEquilibria, ContourStyle, EquilibriumPointColor, Nullclines, NullclineShadingMethod, NullclineStyle} /. {opts} /. Options[NullclinePlot]; contourData = (First /@ Cases[Graphics[ ContourPlot[#, {x, xmin, xmax}, {y, ymin, ymax}, Contours->{0}, DisplayFunction->Identity, Evaluate[FilterOptions[ContourPlot, opts]], PlotPoints->cpp, ContourShading->False]], _Line, Infinity]) & /@ funcs; If[stabQ || (shadedQ && shadmethod === Automatic), stabPts = Check[Equilibria[funcs, {x,xmin, xmax}, {y, ymin, ymax}, contourData[[1]], FilterOptions[Equilibria, opts], EquilibriumFindingMethod -> If[And @@ (FreeQ[funcs, #] & /@ {Abs, Sign, If, Which}), FindRoot, FindMinMethod]], eqFailed]]; Show[If[shadedQ, If[stabPts === eqFailed || shadmethod === SumOfSigns, DensityPlot[{2,1}.Sign[funcs], {x, xmin, xmax}, {y, ymin, ymax}, DisplayFunction->Identity, Mesh->False, PlotPoints-> cpp, PlotRange->{-3, 3}, ColorFunction->(GrayLevel[h[N[#]]]&)], Graphics[{data3 = N@{{{xmin, ymin}, {xmax, ymin}}}; criticalXVals = Join[ If[stabPts == {}, {}, First[Thread[stabPts]]], Union[TurnPts[jj = Join[contourData[[1]], contourData[[2]], data3]]], Flatten[#[[{1,-1}, 1]] & /@ jj], N[{xmin, xmax}] ]; criticalXVals = N[elimDups1D[criticalXVals, 10^-2.5]]; segs = BreakAtXVal[Apply[Join, ToFunctions /@ (elimDups2DOrdered[#, 10^-4] & /@ jj)], criticalXVals]; Table[vv = Select[Select[segs, Length[elimDups1D[ First /@ # , 10^-2.5]] != 1 & ], (* eliminates near-vertical segments *) Abs[#[[1, 1]] - criticalXVals[[i]]] <= 10^-2.5 &]; xxx = (First @ Thread[#] &) /@ vv; mmm = Append[MapIndexed[ (vv[[First[#2], #1]] + vv[[First[#2], #1+1]] )/2 & , (home1[#, (criticalXVals[[i]] + criticalXVals[[i+1]])/2] & /@ xxx)], {(criticalXVals[[i]] + criticalXVals[[i+1]])/2, N@ymax}]; (* here we are matching the segments with their evaluation value, for ease and correctness of sorting *) vvv = Drop[First /@ Sort[Thread[ {Append[vv, {}], mmm}], #1[[-1, -1]] < #2[[-1, -1]] &], -1]; mmm = Sort[mmm, #1[[-1]] < #2[[-1]] &]; grays = GrayLevel /@ (( (Sign[funcs] /. Thread[{x,y}->N[#]]) /. { {1,1}->.2, {0, 0}->1, {1, -1}->.6, {-1,-1}->.4, {-1,1}->.8, {0, _} :> 1, {_, 0} :> 1 }) & /@ ((Drop[mmm, -1] + Rest[mmm])/2) ); Thread[{grays, Polygon[Join[#, {{#[[-1,1]], ymax}, {#[[1,1]], ymax}}]]& /@ vvv}], {i, 1, Length[criticalXVals]-1}]}]], {} ], Graphics[{If[showLines, i = 1; {cs /. Automatic->{}, Append[nullcSty[[1]], Line /@ contourData[[1]]], Append[nullcSty[[2]], Line /@ contourData[[2]]]}, {}], If[stabQ, {stabCol, PointSize[0.021], Point /@ stabPts}, {}] }], Frame->True, FilterOptions[Graphics, opts]]]
If your output will be a sequence of graphs, select a display option.
Display as:
individual images
a movie
both of the above
Up to
Initialization code: The NullclinePlot package