# Basins of Attraction for the Duffing equation with no forcing term

I want to solve an ode for many different initial values and plot these initial values in a plane in

a) red if $(x,\dot{x})$ tends to $(1,0)$ for $t\to\infty$

b) blue if $(x,\dot{x})$ tends to $(-1,0)$ for $t\to\infty$

thanks to the answers to another question I asked, I was able to get this

duffing = x''[t] + 0.02 x'[t] - .5 x[t] + .5 x[t]^3 == 0;
parametricsol =
ParametricNDSolveValue[{duffing, x == x0, x' == v0}, {x, x'}, {t, 0, 4000}, {x0, v0}]

stepsize:=.05

value := Table[
parametricsol[i, j][], {i, -1.5, 1.5, stepsize}, {j, -1.5, 1.5, stepsize}]
position := Table[{i, j}, {i, -1.5, 1.5, stepsize}, {j, -1.5, 1.5, stepsize}]

data := Flatten /@ Transpose@{Catenate@position, Catenate@value}

ListPlot[Style[{#[], #[]},
Which[#[] == -1., White, #[] == 1., Blue, #[] == 0.,
Red]] & /@ data]


If I'm not mistaken ParametricNDSolveValue solves the ode for around 3600 pairs of initial conditions and this takes around 2 minutes on my pc. For a nice looking plot the stepsize=.05 is way too big but if i choose a smaller step size, say .01, it would take almost 1 hour till the computation is done. For a step size of 0.001 the computation would take more than 4 days… What can I do to make this faster?

• Have you seen this? – J. M. is in limbo Jan 29 '17 at 11:35
• thx, no i haven't. Do you know whether there is a discussion of the basin of attraction and how to get it in a fast/nice way in the book? – freddy90 Jan 29 '17 at 12:16
• It's a hint to look in the book (in Google Books or maybe a hard copy borrowed from the nearby library) and find it out for yourself... – J. M. is in limbo Jan 29 '17 at 12:42

# Update

Inspired by the comment by J.M., I managed to code a much faster implementation but to be honest I don't understand the math behind very well:

duffing = x''[t] + 0.02 x'[t] - .5 x[t] + .5 x[t]^3 == 0;
tbegin = -100;
pa2 = ParametricNDSolveValue[{duffing, x == x0, x' == v0},
x, {t, tbegin, 0}, {x0, v0}];

{test, test} = pa2[0, #] & /@ {10^-6, -10^-6};
dat = Transpose[#["ValuesOnGrid"] & /@ {test[#], test[#]'}] & /@ {1, 2};
(* Alternative method for obtaining dat: *)
(*
plot = ParametricPlot[{test[#][t], test[#]'[t]} & /@ {1, 2} // Evaluate, {t, tbegin, 0}]
dat = Cases[plot, Line[a_] :> a, Infinity];
*)

poly = Polygon@
Join[dat[[#2]],
Reverse@dat[[#1]][[Position[dat[[#1]], First@Nearest[dat[[#1]], dat[[#2, 1]]]][[1,
1]] ;;]]] &;

{Red, poly[1, 2], Blue, poly[2, 1], Yellow, Point[{{1, 0}, {-1, 0}}]} // Graphics Add a proper "StopIntegration" event to ParametricNDSolveValue will significantly speed up the calculation:

duffing = x''[t] + 0.02 x'[t] - .5 x[t] + .5 x[t]^3 == 0;
parametricsol =
With[{method = "DetectionMethod" -> "Sign"},
ParametricNDSolveValue[{duffing, x == x0, x' == v0, a == 0,
WhenEvent[x[t] == #, a[t] -> a[t] + #, method] & /@ {-1, 1},
WhenEvent[a[t] == #, Sow@Sign@a[t]; "StopIntegration", method] & /@ {3, -3}},
x, {t, 0, 4000}, {x0, v0}, DiscreteVariables -> a]]

value[{x_: 0}] := x
data = With[{stepsize = 1/100},
ParallelTable[
parametricsol[i, j] // Reap // Last // Flatten // value, {i, -1.5, 1.5,
stepsize}, {j, -1.5, 1.5, stepsize}]]; // AbsoluteTiming
(* {545.779733, Null} dual core *)

ContourPlot[False, {x, -1.5, 1.5}, {xp, -1.5, 1.5},
FrameLabel -> (Style[TraditionalForm@#, 15] & /@ {x, x'})]~Show~
ArrayPlot[data\[Transpose], ColorRules -> {1 -> Red, -1 -> Blue, 0 -> Yellow},
DataRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, DataReversed -> True] • This colors the initial value blue (red) if x is round about -1 (1), right? Should have probably mentioned it but x' has to be round about 0 at the same moment. So I guess StopIntegration doesn't work then? – freddy90 Jan 29 '17 at 13:10
• @freddy90 Managed to make the result more reliable, have a look. – xzczd Jan 29 '17 at 15:06
• Very nice! I also tried to get the separatrix but couldn't figure out how to get the color between the two curves. Thank you! – freddy90 Jan 30 '17 at 11:37
• Great answer. You don't need to use Cases on the graphics object to get the coordinates for each solution. Consider using the methods for InterpolationFunctions instead, like test["ValuesOnGrid"]. – Pillsy Mar 24 '17 at 14:42
• Slightly shorter: dat = Table[Transpose[Through[{test[k], test[k]'}["ValuesOnGrid"]]], {k, 2}]; Also: poly = Polygon[Join[dat[[#2]], dat[[#1]][[-1 ;; First[Nearest[dat[[#1]] -> Automatic, dat[[#2, 1]]]] ;; -1]]]] &; – J. M. is in limbo Mar 25 '17 at 8:24

Comment

I fail to see the point of taking a very very small step size. Even, with a few points you can get a reasonable plot.

sol[x0_?NumericQ]:=First@NDSolve[{x''[t] + 0.02 x'[t] - .5 x[t] + .5 x[t]^3 == 0,
x == x0, x' == x0}, {x}, {t, 0, 4000}];
With[{sol1 = {x[t], x'[t]} /. sol[#] & /@ Range[-1.5, 1.5, 0.05]},
ParametricPlot[sol1, {t, 0, 4000}, PlotRange -> All,PlotPoints -> 500,
MaxRecursion -> 5,AxesLabel -> {"x", "x'"}]] // AbsoluteTiming The next thing come to mind was to convert your second order ode to a system of first order odes, which happens to be a little faster than the first one.

soln[x0_?NumericQ] :=First@NDSolve[{x1'[t] == x2[t],
x2'[t] + 0.02 x2[t] - .5 x1[t] + .5 x1[t]^3 == 0, x1 == x0,
x2 == x0}, {x1, x2}, {t, 0, 4000}];
With[{sol2 = {x1[t], x2[t]} /. soln[#] & /@ Range[-1.5, 1.5, 0.05]},
ParametricPlot[sol2, {t, 0, 4000}, PlotRange -> All,
PlotPoints -> 500, MaxRecursion -> 5,
AxesLabel -> {"x1", "x2"}]] // AbsoluteTiming A much faster way to do it, is to plot the stream lines and then impose the solution curves,

sp = StreamPlot[{x2, -0.02*x2 + 0.5*x1 - 0.5*x1^3}, {x1, -2,
2}, {x2, -2, 2}, StreamColorFunction -> "Rainbow"];
Manipulate[
Show[sp,
ParametricPlot[
Evaluate[
First[{x1[t], x2[t]} /.
NDSolve[{x1'[t] == x2[t],
x2'[t] + 0.02 x2[t] - .5 x1[t] + .5 x1[t]^3 == 0,
Thread[{x1, x2} == point]}, {x1, x2}, {t, 0, T}]]], {t,
0, T}, PlotStyle -> Red]], {{T, 4000}, 1,
100}, {{point, {-1.5, 1.5}}, Locator}, SaveDefinitions -> True] I am thankful to @halirutan for teaching me how to use AbsoluteTiming.

• I'm afraid this answer doesn't meet the requirement of OP, he wants a) red if the solution tends to $1$ for $t→∞$ b) blue if the solution tends to $-1$ for $t→∞$ – xzczd Jan 29 '17 at 8:22
• @xzczd hmm... You are right. This is not a complete answer. My focus was on "What can I do to make this faster?" – zhk Jan 29 '17 at 8:26