# Simple way to highlight streams in basins of attraction in StreamDensityPlot

In making a figure to answer this question, I wanted to highlight streams that start or end at the critical points, $$(0,0)$$ and $$(2,2)$$.

One can define StreamPoints or VectorPoints, but that doesn't create the streams both to and away from a critical point. The only other way seems to be rather awkward: making a ParametricPlot and superimposing it on the StreamDensityPlot.

Question: How can I most simply alter the below code to show (in red) all streams originating from the local optimum at $$(2,2)$$ and (in green) all streams leaving or terminating at the saddle point at $$(0,0)$$? (Some stream lines will be "both"... i.e., leave $$(2,2)$$ and terminate at $$(0,0)$$.)

StreamDensityPlot[{{3 x^2 - 6 y, 3 y^2 - 6 x}, x^3 + y^3 - 6 x y},
{x, -5, 5}, {y, -5, 5},
StreamPoints -> {{{{1, 0}, Red}, {{-1, -1}, Green}, Automatic}},
Epilog -> {Red, PointSize[0.03], Point[{2, 2}], Green,
Point[{0, 0}]}]


This is a partial solution:

myStreams =
Table[{{2, 2} + 2 {Cos[θ], Sin[θ]}, Red}, {θ, 0, 2 π, .3}];
StreamDensityPlot[{{3 x^2 - 6 y, 3 y^2 - 6 x},
x^3 + y^3 - 6 x y},
{x, -5, 5}, {y, -5, 5},
StreamPoints ->
Flatten[{Join[myStreams, {{{-.2, -.2}, Green}}], Automatic}, 1],
Epilog -> {Red, PointSize[0.03], Point[{2, 2}], Green, Point[{0, 0}]}]


• For the green streams going in and out of the saddle point, do you just want the stable and unstable manifolds (4 lines -- in from 2 sides and out from 2 others)? Commented Dec 30, 2018 at 23:40
• Four lines would suffice. Commented Dec 31, 2018 at 0:08

We could get the positions of the tips of the arrows from the graphics itself, integrate to see where test points at those positions end up, and then color the arrows accordingly. Here is the code for that:

plot = StreamDensityPlot[
{{3 x^2 - 6 y, 3 y^2 - 6 x}, x^3 + y^3 - 6 x y},
{x, -5, 5}, {y, -5, 5},
Epilog -> {Red, PointSize[0.03], Point[{2, 2}], Green, Point[{0, 0}]}
]

arrows = Cases[plot, _Arrow, Infinity];
tips = arrows[[All, 1, -1]];

headedIntoBasin = intoBasinQ[{0, 0}, #] & /@ tips;
headedFromBasin = fromBasinQ[{2, 2}, #] & /@ tips;

plot /. Join[
# -> {Yellow, #} & /@ Pick[arrows, both],
# -> {Green, #} & /@ Pick[arrows, headedIntoBasin],
# -> {Red, #} & /@ Pick[arrows, headedFromBasin]
]


The functions intoBasinQ and fromBasinQ are verbose so I leave them for last, although they are quite simple, they only look complicated:

intoBasinQ[basin_, {x0_, y0_}] := Module[{xfun, yfun},
{xfun, yfun} = Quiet@NDSolveValue[{
x'[t] == 3 x[t]^2 - 6 y[t],
y'[t] == 3 y[t]^2 - 6 x[t],
x[0] == x0,
y[0] == y0,
WhenEvent[
Norm[basin - {x[t], y[t]}] < 0.1,
"StopIntegration",
"LocationMethod" -> "StepEnd"
]
}, {x, y}, {t, 0, 10}];
xf = Last@Flatten@xfun["ValuesOnGrid"];
yf = Last@Flatten@yfun["ValuesOnGrid"];
Norm[basin - {xf, yf}] < 0.2
]

fromBasinQ[basin_, {x0_, y0_}] := Module[{xfun, yfun},
{xfun, yfun} = Quiet@NDSolveValue[{
x'[t] == 3 x[t]^2 - 6 y[t],
y'[t] == 3 y[t]^2 - 6 x[t],
x[0] == x0,
y[0] == y0,
WhenEvent[
Norm[basin - {x[t], y[t]}] < 0.1,
"StopIntegration",
"LocationMethod" -> "StepEnd"
]
}, {x, y}, {t, 0, -10}];
xf = First@Flatten@xfun["ValuesOnGrid"];
yf = First@Flatten@yfun["ValuesOnGrid"];
Norm[basin - {xf, yf}] < 0.2
]

• Oh how nice. Just what I needed. (accept) Perhaps Wolfram will include some of this functionality. Commented Dec 31, 2018 at 1:08

Similar idea to @C.E.'s, but using StreamColorFunction, which flummoxed me, since it does not work as documented for StreamDensityPlot, when the argument is of the form {vector field, scalar field}:

vf2ode[vf_, vars_List] :=  (* vector field to ode *)
D[Through[vars@t], t] == (vf /. Thread[vars -> Through[vars@t]]);

(* StreamColorFunction *)
myColor[0. | 0] = ColorData[97][1];
myColor[1. | 1] = Green;
myColor[2. | 2] = Red;
myColor[3. | 3] = Purple;  (* hits both singular points*)
myColor[_] = Black;        (* shouldn't happen *)

scf = Function[{xx, yy},  (* stream color function *)
Which[
Norm[{xx, yy} - {0., 0.}] < 10^-8, myColor[1.],
Norm[{xx, yy} - {2., 2.}] < 10^-8, myColor[2.],
True, myColor@Total[
Block[{x, y, t, color},
NDSolveValue[{
vf2ode[{3 x^2 - 6 y, 3 y^2 - 6 x}, {x, y}], {x[0], y[0]} == {xx, yy},
color[0] == 0,
WhenEvent[Abs[x[t]] > 5.1, "StopIntegration"],
WhenEvent[Abs[y[t]] > 5.1, "StopIntegration"],
WhenEvent[Norm[{x[t], y[t]} - cp[[1]]] < 10^-1, (* unstable => large tol. *)
{color[t] -> color[t] + 1, "StopIntegration"}],
WhenEvent[Norm[{x[t], y[t]} - cp[[2]]] < 10^-4,
{color[t] -> color[t] + 2, "StopIntegration"}]},
color["ValuesOnGrid"],
{t, -100, 100},
StartingStepSize -> 0.001,
DiscreteVariables -> {color}
][[{1, -1}]]
]
]
]
];

(* unstable separatrices *)
sp = Map[Last,
NDSolveValue[{
vf2ode[{3 x^2 - 6 y, 3 y^2 - 6 x}, {x, y}], {x[0], y[0]} == #,
WhenEvent[Abs[x[t]] > 3.3, "StopIntegration"],
WhenEvent[Abs[y[t]] > 3.3, "StopIntegration"]},
{x["ValuesOnGrid"], y["ValuesOnGrid"]},
{t, 0, 100},
StartingStepSize -> 0.001, PrecisionGoal -> 10, AccuracyGoal -> 15
] & /@ ({{-1, 1}, {1, -1}}/10^8),
{2}]


Graphics:

Show[
DensityPlot[x^3 + y^3 - 6 x y,
{x, -5, 5}, {y, -5, 5},
Epilog -> {Red, PointSize[0.03], Point[{2, 2}], Green, Point[{0, 0}]},
PlotRange -> All],
StreamPlot[{3 x^2 - 6 y, 3 y^2 - 6 x},
{x, -5, 5}, {y, -5, 5},
StreamPoints -> {{{1, 1}, {3, 3}, {-1, -1}, Sequence @@ sp, Automatic}},
StreamColorFunction -> scf, StreamColorFunctionScaling -> False]
]


Just for fun here's the lift of the mapping $$S^2 \rightarrow {\Bbb{RP}^2}$$ of the phase portrait on the real projective plane of the projectivization of the ODE. Antipodal points of the sphere $$S^2$$ should be identified to get $${\Bbb{RP}^2}$$. It can be projectivized because the vector field, which is polynomial, can easily be made homogeneous. We can see there's another critical point at infinity $$[x \colon y \colon z] = [1 \colon -1 \colon 0]$$. This c.p. becomes more apparent in the StreamPlot if the domain is extended to {x, -100, 100}, {y, -100, 100}: The slopes of the stream lines where they intersect $$y = -x$$ approach 1 at infinity. Code here: https://pastebin.com/84dTTbHs

• Superb. Thanks so much. (+1) Wolfram should include this functionality. Commented Dec 31, 2018 at 1:46