8
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ 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)? $\endgroup$ – Chris K Dec 30 '18 at 23:40
  • $\begingroup$ Four lines would suffice. $\endgroup$ – David G. Stork Dec 31 '18 at 0:08
14
$\begingroup$

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

Mathematica graphics

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
  ]
$\endgroup$
  • 1
    $\begingroup$ Oh how nice. Just what I needed. (accept) Perhaps Wolfram will include some of this functionality. $\endgroup$ – David G. Stork Dec 31 '18 at 1:08
10
$\begingroup$

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

enter image description here


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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Superb. Thanks so much. (+1) Wolfram should include this functionality. $\endgroup$ – David G. Stork Dec 31 '18 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.