Basins of Attraction

Question

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.

Answer 1

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

Answer 2

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