# StreamPoints for all of the domain (how to color basins of attraction)

Although that title is probably as confusing as I'll get out, here is the actual question. In StreamPlot, the StreamPoints will help trace out a 'trajectory'. So this nice, but I'm interested in areas throughout the domain that will converge to a certain location.

So for the following function, (often call the unforced Duffing system) I would like to have a StreamPlot and in the background almost like a contour plot of two colors (I know a contour plot isn't right at all I'm just trying to help you see what I'm interested in). So if it is one color then it will converge at one hole and another color for the other hole. Not sure if this is possible, but any ideas would be awesome!

StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
StreamPoints -> {{{{1, 1}, Red}, {{1, 1.6}, Green}, Automatic}}]


• maybe StreamDensityPlot, e.g. StreamDensityPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}, StreamPoints -> {{{{1, 1}, Red}, {{1, 1.6}, Green}, Automatic}}, ColorFunction -> "Rainbow"]? – kglr Feb 16 '18 at 22:57
• @kglr Sorry, I thought it might be a bad idea to include contour plot at all. So the StreamDensityPlot just shows the strength of the StreamPlot (much like a VectorPlot would just have a larger vector). I basically want a background (red and green for example), and at say (1,1) it would be green if the Streamline was going to converge to the right hole. And I want that for the whole domain. Once again, this might not be possible. – Josh Feb 16 '18 at 23:20
• Sounds like these q/a s may be useful. – kglr Feb 16 '18 at 23:58
• Closely related: Fractal basins of attraction in a Magnetic Pendulum. – C. E. Feb 26 '18 at 5:05

Horribly slow, but you probably get the idea. If this system can be solved analytically, then one could implement a far better approach, but I haven't looked into this.

What you are after is to find out which point in your plane converges to which sink. You can do this by solving the differential equation where your direction vector is the speed. You provide a starting point and start integrating. If you have done this, you can look what point is reached at the end of the integration and use the specific color for the sink.

A function that does exactly this is the following. It additionally stores the path starting from the initial point to the end.

line[{x0_, y0_}] :=
Module[{x, y, t, res, col},
res = NDSolveValue[{x'[t] == y[t],
y'[t] == -.15 y[t] + x[t] - x[t]^3, x[0] == x0, y[0] == y0}, {x,
y}, {t, 0, 100}];
col = If[Norm[{-1, 0} - Through[res[100]]] < .01, Green, Red];
{col, Line[Table[Through[res[t]], {t, 0, 100, .1}]]}
]


Now you can sample your region and for each starting point you create a colored line (this may take a while)

lines = ParallelTable[line[{x, y}], {x, -2, 2, .1}, {y, -2, 2, .1}];


After that, you can combine all lines into one graphics. This is again very slow due to Opacity and the sheer amount of line points

Show[
Graphics[{Opacity[.01], Thickness[0.01], lines}],
StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}],
Frame -> True, FrameTicks -> True, PlotRange -> {{-2, 2}, {-2, 2}},
PlotRangeClipping -> True
]


but the result looks somewhat nice and I believe it is what you are after

Please be aware that there is much room for improvement.

## Line integral convolution

Another idea is to simply provide a test-function that tells you to which sink the starting point converges. This is close to the line function but doesn't create so many line-points.

test[{x0_?NumericQ, y0_?NumericQ}] :=
Module[{x, y, t, res},
res = NDSolveValue[{x'[t] == y[t],
y'[t] == -.15 y[t] + x[t] - x[t]^3, x[0] == x0, y[0] == y0}, {x,
y}, {t, 0, 100}];
Norm[{-1, 0} - Through[res[100]]] < .01
]


Now, you could create a colored image that shows you the regions of convergence in different colors. We will split a color scheme into two halves and introduce some noise which will be an excellent starting point for a LIC

mat = ParallelTable[
List @@
ColorData["TemperatureMap",
Boole[test[{x, y}]]/2 + RandomReal[1/2]], {y, -2,
2, .01}, {x, -2, 2, .01}];

bg = LineIntegralConvolutionPlot[{{y, -.15 y + x - x^3},
Image[Reverse@mat]}, {x, -2, 2}, {y, -2, 2}, RasterSize -> 512]


and then

Show[
bg,
StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
StreamStyle -> Gray]
]


and finally, note that you can use this test-function with RegionPlot or you can employ it to color your StreamPlot accordingly:

stP = StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
StreamColorFunction ->
Function[{x, y},
If[test[{x, y}], ColorData["TemperatureMap", .9],
ColorData["TemperatureMap", .1]]],
StreamColorFunctionScaling -> False]

Show[bg, stP]


## Speed

As mentioned in the comments, the real system has no analytical solution. The bottlenecks in the last LIC example are (1) the creation of the domain image and (2) the calculation of the LIC itself. As for (2) there is nothing much we can do. However, one can speed up (1) by implementing a simple numerical scheme that can be compiled. Here is a hack of a simple Runge-Kutta solver that can be called in parallel. If you don't have a C compiler then you need to remove the compile to C option:

f[{x_, y_}] := {y, -.15 y + x - x^3}
step[yn_, h_] := Module[{k1, k2, k3},
k1 = f[yn];
k2 = f[yn + h/2.0*k1];
k3 = f[yn - h*k1 + 2.0 h*k2];
yn + h (1./6 k1 + 4./6 k2 + 1./6 k3)
]

fc = Compile[{{x0, _Real, 0}, {y0, _Real, 0}},
Module[{eps = 0.001, xn = x0, yn = y0, h = .1, maxSteps = 5000,
sink = -1, xtmp = 0., ytmp = 0.},
Do[
{xtmp, ytmp} = #;
If[
Norm[{xtmp, ytmp} - {xn, yn}] < eps,
Break[]
];
{xn, yn} = {xtmp, ytmp},
maxSteps
];
If[Norm[{-1, 0} - {xn, yn}] < .1, 1, 0]
], Parallelization -> True, RuntimeAttributes -> {Listable},
CompilationTarget -> "C"
] &[FullSimplify[step[{xn, yn}, h]]]


For a region sampling of 512x512 points, the runtime of the call to NDSolve within the test function is about 100s, the same compution with fc only needs 4s.

mat = fc @@ Transpose[
Table[{x, y}, {y, -2, 2, 4/511}, {x, -2, 2, 4/511}], {2, 3,
1}]; // AbsoluteTiming
Image[Reverse[mat]]


I promise it will be the last graphics, but I quite like the style. It's in remembrance of a good friend of mine, who had a fast line integral convolution for Mathematica long before Wolfram could even spell the word. He always liked flow fishes in vector fields (which are called "Drop" in Mathematica) and often combined them:

lic = LineIntegralConvolutionPlot[{{y, -.15 y + x - x^3},
GaussianFilter[Reverse[mat + RandomReal[{-1, 1}, {512, 512}]],
4]]]},  {x, -2, 2}, {y, -2, 2}, RasterSize -> 512];
Show[
lic,
StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
StreamStyle -> White, StreamMarkers -> "Drop"]
]


• Haha, this is pretty good and it does do the job. I was hoping there could be a way around doing this though! I did mark it as useful, but I'm holding off for a bit to accept it to see if anyone else comes up with something. Thanks! – Josh Feb 17 '18 at 2:25
• @Josh I added certain additional ideas. Would you know if there is a closed form of the differential equation? – halirutan Feb 17 '18 at 2:55
• I just did the first method above, but changed it up a bit. I ran the integration backwards, changed the boolean statement, and then just ran the 'line' function just through the center of the sink (that way you don't have as many lines you have to plot), – Josh Feb 17 '18 at 2:56
• Actually, this isn't the equation that I am working on (the one I'm working on is nonlinear and 4 equations instead of two, and is definitely not closed form). I just wanted to see the concept on something simple. Thanks though. – Josh Feb 17 '18 at 3:01
• @Josh Ahh, good. At least someone who is capable of providing a good minimal example. It seems rare these days. In that case, your goal should be to find a small and fast integrator for your equations that has enough accuracy to converge to the correct sink and can probably be Compiled and ran in parallel. – halirutan Feb 17 '18 at 3:09

In Jumps in NDSolve results, I found the separatrices and used them to construct polygons to color the basins of attraction, but since it had little to do with the question asked, I left the code for it out. Below is a similar approach. The plot is just two polygons, clipped to the plot range.

eq = {y, -15/100 y + x - x^3} /. {v : x | y :> v[t]}  (* ode system *)
(*  {y[t], x[t] - x[t]^3 - (3 y[t])/20}  *)

equil = Solve[eq == 0]   (* equilibria *)
(*  {{x[t] -> -1, y[t] -> 0}, {x[t] -> 0, y[t] -> 0}, {x[t] -> 1, y[t] -> 0}}  *)

(* Picking the separatrices that are the boundaries of the basins takes some
* inspection of the system: Reverse the velocity field (eq -> -eq), and the basins
* should be the same; however, the eigenvalues below would change sign.
*)
sepeq = Pick[N@equil, Det@D[eq, {{x[t], y[t]}}] /. equil // Sign, -1] (* unstable ones *)
sep0 = Flatten[  (* initial velocity vectors of separatrices *)
{#, -#} &@ # & /@
Pick[#2, Sign[#1], -1] &,  (* picks eigenvectors heading toward equilibrium *)
Eigensystem[#]
],
1] & /@ (D[eq, {{x[t], y[t]}}] /. sepeq)
(*
{{x[t] -> 0., y[t] -> 0.}}
{{{-0.680151, 0.680151}, {0.733072, -0.733072}}}
*)

ics = MapThread[ (* IC = initial velocity times a small number *)
{x[t], y[t]} + 100 # 2 $MachineEpsilon & /@ #2 /. #1 &, {sepeq, sep0} ]; sepsols = Flatten[ (* separatrix integral curves *) Map[ NDSolve[{{x'[t], y'[t]} == eq, {x[0], y[0]} == #, WhenEvent[Abs[x[t]] > 2.2 && Abs[y[t]] > 2.2, "StopIntegration"]}, {x, y}, {t, 0, -100}, (* eigenvalue is neg. so integrate backwards *) "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}] &, ics, {2} ], 2]; lines = Cases[ (* get points of separatrices *) Normal@ ParametricPlot[Evaluate[{x[t], y[t]} /. sepsols], {t, -100, 100}, PlotPoints -> 100, PlotRange -> All], Line[p_] :> p, Infinity]; basins = Module[{i}, (* boundaries of basins of attraction *) i = First@Nearest[#[[2]] -> "Index", Last@#[[1]]]; Join[#[[1]], #[[2, i ;; 1 ;; -1]]] ] & /@ {Reverse /@ lines, Reverse[Reverse /@ lines]}; Show[ Graphics[{Opacity[0.5], Riffle[{Red, Green}, Polygon /@ basins]}], #, Options@# ] &@ StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2}, StreamPoints -> {{{{1, 1}, Darker@Red}, {{1, 1.6}, Darker@Green}, Automatic}}, PlotRange -> 2.05]  Update: We can use the other separatrix as one of the stream lines and add a few more colored appropriately. sep2 = Flatten[{#, -#} &@# & /@ MapThread[Pick[#2, Sign[#1], 1] &, Eigensystem[#]], 1] & /@ (D[eq, {{x[t], y[t]}}] /. sepeq); ics2 = MapThread[{x[t], y[t]} + 100 # 2$MachineEpsilon & /@ #2 /. #1 &, {sepeq, sep2}];
sepsols2 = Flatten[
Map[
NDSolve[{{x'[t], y'[t]} == eq, {x[0], y[0]} == #,
WhenEvent[x[t]^2 + y[t]^2 > 0.25, "StopIntegration"]},
{x, y}, {t, 0, 100}] &,
ics2, {2}],
2];

Show[Graphics[{Opacity[0.5],
Riffle[{Red, Green}, Polygon /@ basins]}], #, Options@#] &@
StreamPlot[{y, -.15 y + x - x^3}, {x, -2, 2}, {y, -2, 2},
StreamPoints -> {Join[