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 can I plot the stable and unstable separatrices to make a graphic of the basins of attraction of this differential equation?

y'[t] == Cos[t] - Sin[x[t]] - 0.1 y[t], x'[t] == y[t]

I used as a reference this paper from hubbard and I think I can draw the different basins of attraction using separatrices. I'm using this Poincarè map:

P[{a_, b_}, tmax_] := {x[tmax], y[tmax]} /. First[NDSolve[
Join[{y'[t] == Cos[t] - Sin[x[t]] - 0.1 y[t], 
  x'[t] == y[t]}, {x[0] == a, y[0] == b}], {x, y}, {t, 0, tmax}]]

and I did plot one of the basins as points with this:

(*It will take a little bit of time to see the graphic*)
test = NestList[P[#, 4 \[Pi]] &, {1, 1}, 5000];
cb = Last[test]
pointlist = 
  Compile[{{a, _Real}, {b, _Real}}, 
   FixedPointList[P[#, 2 \[Pi]] &, {a, b}, 8000, 
    SameTest -> (Abs[#1[[1]] - #2[[1]]] < 0.1 && 
        Abs[#1[[2]] - #2[[2]]] < 0.1 &)]];
listofpoint = 
   Cases[Flatten[
     Table[If[
       Abs[cb[[1]] - Last[pointlist[a, b]][[1]]] < 0.1 && 
        Abs[cb[[2]] - Last[pointlist[a, b]][[2]]] < 0.1, 
       pointlist[a, b], {x}], {a, -5, 5, 0.1}, {b, -5, 5, 0.1}], 
     1], Except[{x}]];
ListPlot[listofpoint, PlotStyle -> {PointSize[Tiny], Red}, 
 Background -> Black]

I think that I have to find the saddle points to draw at least one basin but I'm stuck. On a side note I used the following method for the Duffing equation but in this case it didn't work:

Duffing = {x'[t] == y[t], y'[t] == 0.5 x[t] - 0.5 x[t]^3 - 0.015 y[t]};
{x, y} /. 
 Solve[Thread[(Last /@ Duffing /. a_[t] :> a) == {0, 0}], {x, y}]
separatrix = ({x[t], y[t]} /. 
      First[NDSolve[{Duffing, x[0] == #[[1]], y[0] == #[[2]]}, {x[t], 
         y[t]}, {t, -133, 0}]] &) /@ {{0, -0.001}, {0, 0.001}};
sepPlot = 
 Show[ParametricPlot[Evaluate[separatrix[[1]]], {t, -132.6, 0}], 
  ParametricPlot[Evaluate[separatrix[[2]]], {t, -129.5, 0}], 
  Axes -> None, Frame -> True]
data = Cases[sepPlot, _Line, \[Infinity]] /. Line -> List;
Short[data, 5]
Graphics[{Yellow, Rectangle[{-1.7, -1}, {1.7, 1}], Blue, 
  Polygon[Join[data[[1, 1]], Reverse[data[[2, 1]]]]], 
  PointSize[0.02], {GrayLevel[1], Point[{1, 0}]}, Point[{-1, 0}]}, 
 Frame -> True, FrameTicks -> {Automatic, Automatic, None, None}, 
 PlotRange -> {{-1.7, 1.7}, {-1, 1}}, PlotRangeClipping -> True]

P.S:I'm using Mathematica 8.0

share|improve this question
    
Is this question somehow related to your question? –  Rod Aug 26 '13 at 3:16
    
@rod-lm yes it's related. The first answer that uses RegionPlot could solve my problem, but I tried it for my equation and it didn't work –  Iulian Aug 27 '13 at 13:31
add comment

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.