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
RegionPlotcould solve my problem, but I tried it for my equation and it didn't work – Iulian Aug 27 '13 at 13:31