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$
  • 2
    $\begingroup$ Have you seen this? $\endgroup$ – J. M. will be back soon Jan 29 '17 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 '17 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$ – J. M. will be back soon Jan 29 '17 at 12:42
9
$\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$
  • $\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 '17 at 13:10
  • 1
    $\begingroup$ @freddy90 Managed to make the result more reliable, have a look. $\endgroup$ – xzczd Jan 29 '17 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 '17 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 '17 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$ – J. M. will be back soon Mar 25 '17 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
    $\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 '17 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 '17 at 8:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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