5
$\begingroup$

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[0] == x0, x'[0] == v0}, {x, x'}, {t, 0, 4000}, {x0, v0}]

stepsize:=.05    

value := Table[
parametricsol[i, j][[1]][3000], {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[{#[[1]], #[[2]]}, 
Which[#[[3]] == -1., White, #[[3]] == 1., Blue, #[[3]] == 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?

$\endgroup$
3
  • 2
    $\begingroup$ Have you seen this? $\endgroup$ Jan 29, 2017 at 11:35
  • $\begingroup$ 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? $\endgroup$
    – freddy90
    Jan 29, 2017 at 12:16
  • 1
    $\begingroup$ 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... $\endgroup$ Jan 29, 2017 at 12:42

2 Answers 2

10
$\begingroup$

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[0] == x0, x'[0] == v0}, 
   x, {t, tbegin, 0}, {x0, v0}];

{test[1], test[2]} = 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

Mathematica 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[0] == x0, x'[0] == v0, a[0] == 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]

Mathematica graphics

$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – freddy90
    Jan 29, 2017 at 13:10
  • 1
    $\begingroup$ @freddy90 Managed to make the result more reliable, have a look. $\endgroup$
    – xzczd
    Jan 29, 2017 at 15:06
  • $\begingroup$ 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! $\endgroup$
    – freddy90
    Jan 30, 2017 at 11:37
  • 2
    $\begingroup$ 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[1]["ValuesOnGrid"]. $\endgroup$
    – Pillsy
    Mar 24, 2017 at 14:42
  • 1
    $\begingroup$ 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]]]] &; $\endgroup$ Mar 25, 2017 at 8:24
5
$\begingroup$

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[0] == x0, x'[0] == 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

enter image description here

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[0] == x0, 
     x2[0] == 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

enter image description here

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[0], x2[0]} == point]}, {x1, x2}, {t, 0, T}]]], {t, 
    0, T}, PlotStyle -> Red]], {{T, 4000}, 1, 
  100}, {{point, {-1.5, 1.5}}, Locator}, SaveDefinitions -> True]

enter image description here

I am thankful to @halirutan for teaching me how to use AbsoluteTiming.

$\endgroup$
2
  • 2
    $\begingroup$ 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→∞$ $\endgroup$
    – xzczd
    Jan 29, 2017 at 8:22
  • $\begingroup$ @xzczd hmm... You are right. This is not a complete answer. My focus was on "What can I do to make this faster?" $\endgroup$
    – zhk
    Jan 29, 2017 at 8:26

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