# Basins of Attraction

How does one shade the basin(s) of attraction of a phase plot in Mathematica? I have been trying to do this using the system

\begin{align*} \dot x &= y\\ \dot y &= -9\sin(x) - 0.20y \end{align*}

but I have gotten nowhere.

-

Here's a bruteforce way to do it for the simple case when the attractor is a fixed point.

1. Find a fixed point
2. Pick initial values for ODE
3. Solve ODE, see if it gets close to fixed point
4. Go back to 2 until satisfied

By looking at Reduce[y == 0 && -9 Sin[x] - 2/10 y == 0, {x, y}, Reals] we see that for instance {x=0,y=0} is a fixed point, let's use that

tmax = 40;
tol = 0.2;

(* Solution to ODE that maps t to {x[t],y[t]} *)
sol[x0_?NumericQ, y0_?NumericQ] := First@NDSolve[{
x[0] == x0,
y[0] == y0,
x'[t] == y[t],
y'[t] == -9 Sin[x[t]] - 0.2 y[t]},
{x, y}, {t, 0, tmax}] /.
HoldPattern[{x -> xi_, y -> yi_}] :>
Function[{t}, {xi[t], yi[t]}];

(* Create function that takes a solution as argument and
checks if it's close to attractor at tmax*)
MakeBasinTest[{x0_, y0_}] := Function[{f}, Norm[f[tmax] - {x0, y0}] <= tol];

inCenterBasinQ = MakeBasinTest[{0, 0}]

(* Create streamplot and basin region, give it a while *)
Show[
StreamPlot[{y, -9 Sin[x] - 0.2 y}, {x, -10, 10}, {y, -10, 10}],
RegionPlot[
inCenterBasinQ[sol[x0, y0]], {x0, -10, 10}, {y0, -10, 10},
PlotPoints -> 30, MaxRecursion -> 1,
PlotStyle -> {Green, Opacity[0.2]}]]


You can change tmax and the options for RegionPlot until you have a good time/quality tradeoff

-
Take a look at ParametricNDSolveValue. I used pf = With[{tmax = 100}, ParametricNDSolveValue[{x'[t] == y[t], y'[t] == -9 Sin[x[t]] - 0.20 y[t], x[0] == a, y[0] == b}, {x[tmax], y[tmax]}, {t, 0, tmax}, {a, b}]] to simplify the first part of your code. –  Szabolcs Apr 27 at 18:32
@Szabolcs The ParametricNDSolve stuff is neat, alas I only have v8 :) Was there any significant improvment in speed? –  ssch Apr 27 at 19:49

No "brute-force" playing with NDSolve, we can get an idea of attraction basins with the StreamDensityPlot and StreamPoints option in it.

Let's find e few points of interest where the vector flow becomes zero. E.g.

{x, y} /. {ToRules @ LogicalExpand @
Reduce[y == 0 && -9 Sin[x] - 1/5 y == 0 && -5 < x < 10, {x, y}]}

 {{0, 0}, {-Pi, 0}, {Pi, 0}, , {2 Pi, 0}, {3 Pi, 0}}


We are looking what happens near points close to these solutions i.e. {{-Pi, 0}, {Pi, 0}, {3 Pi, 0}}. Define an epsilon :

eps = 1/20;


Now we can observe with StreamDensityPlot behavior of appropriate points. We could change eps with e.g. Manipulate, here we demonstrate standard output :

GraphicsColumn[
Table[
StreamDensityPlot[{y, -9 Sin[x] - 1/5 y}, {x, -12, 12}, {y, -8, 8},
StreamPoints -> {{
{{-Pi + k eps, eps}, Directive[Thick, Red]},
{{Pi + k eps, eps},  Directive[Thick, Darker @ Green]},
{{3 Pi + k eps, -eps}, Directive[Thick, Orange]}, Automatic}},
AspectRatio -> Automatic, ImageSize -> 500],
{k, {-1, 1}}]]


These plots give quite a good idea of attraction basins

-