Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

2 Answers 2

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

basin

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

share|improve this answer
1  
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 '13 at 18:32
    
@Szabolcs The ParametricNDSolve stuff is neat, alas I only have v8 :) Was there any significant improvment in speed? –  ssch Apr 27 '13 at 19:49
add comment

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

enter image description here

These plots give quite a good idea of attraction basins

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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