Qfyilyi

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}]}]

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;

both = Thread[headedIntoBasin && headedFromBasin];



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

  ]

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: needs 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]

 ]