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.
$\begingroup$ $\endgroup$
1 
- $\begingroup$ How would you shade strange attractors - of which there are various types. $\endgroup$alancalvitti– alancalvitti2014-11-05 18:19:39 +00:00Commented Nov 5, 2014 at 18:19
Add a comment |
2 Answers
$\begingroup$ $\endgroup$
3 Here's a bruteforce way to do it for the simple case when the attractor is a fixed point.
- Find a fixed point
- Pick initial values for ODE
- Solve ODE, see if it gets close to fixed point
- 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
- 1$\begingroup$ Take a look at
ParametricNDSolveValue. I usedpf = 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. $\endgroup$Szabolcs– Szabolcs2013-04-27 18:32:01 +00:00Commented Apr 27, 2013 at 18:32 - $\begingroup$ @Szabolcs The
ParametricNDSolvestuff is neat, alas I only have v8 :) Was there any significant improvment in speed? $\endgroup$ssch– ssch2013-04-27 19:49:48 +00:00Commented Apr 27, 2013 at 19:49 - $\begingroup$ I cannot make it work for
StreamPlot[{y, -0.5 y - Sin[x]}, {x, -3 \[Pi], 3 \[Pi]}, {y, -2, 2}, ....](changed insolas well), withtmax=500andPlotPoints -> 200, MaxRecursion -> 5, after waiting for a long time I got just a hexagonal-like box. Any ideas how to make it more flexible? $\endgroup$corey979– corey9792021-04-11 11:37:17 +00:00Commented Apr 11, 2021 at 11:37
$\begingroup$ $\endgroup$
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
