Sweep-up for Duffing equation

I would like to calculate a sweep-up for the Duffing equation. To be more precise, I'm looking for a solution to the ode

$\qquad x''(t)+0.1x'(t)+x(t) +0.01 x^3(t) = 0.8 \cos(\eta t)$

Let's say for the initial conditions $x(0)=x'(0)=0$ and $\eta \in [0.8, 1.2]$ with increment size $0.001$ and $t \in [0,2\cdot 100 \frac{\pi}{\eta}]$. But instead of using the same initial conditions for every $\eta$, I would like to choose new initial conditions $x_{j+1}(0)=x_j(2\cdot 100 \frac{\pi}{\eta})$ and $x_{j+1}'(0)=x_j'(2\cdot 100 \frac{\pi}{\eta})$ every time $\eta$ increases, where $x_j$ is the solution of the previous computation with the previous $\eta$. As output I would like to have the maximum value of $x$ for every $\eta$ neglecting the transients.

I tried it using Table, which yields reasonable results but takes some time to calculate. Here is the code

numbper = 100; (*number of periods*)
sol = ParametricNDSolveValue[{x''[t] + 1/10 x'[t] + x[t] + 1/100 x[t]^3 ==
8/10 Cos[η t], x[0] == x0, x'[0] == v0}, {x, x'}, {t, 0, numbper*2 Pi/η},
{x0, v0, η}];

xstart = sol[0, 0, 0.79][[1]][numbper*2 Pi/0.79];
vstart = sol[0, 0, 0.79][[2]][numbper*2 Pi/0.79];
dη = 0.001;
sweepup = Table[{
NMaxValue[{sol[xstart, vstart, η][[1]][
t], (numbper - 1)*2 Pi/η <= t <= numbper*2 Pi/η},
t], xstart =
sol[xstart, vstart, η][[1]][(numbper - 1)*2 Pi/η];
vstart = sol[xstart, vstart, η][[2]][numbper*2 Pi/η];}
, {η, 0.8, 1.2, dη}];


I'm now wondering what I could do to speed things up.

• My advice: Better use Do instead of For. Write it down as it comes to your mind and optimize later. Jan 11, 2018 at 17:54
• thanks, i will try that out. that's probably a dumb question but is it even advisable to use FoldList here? I mean, can we save computation time by doing so? Jan 11, 2018 at 18:32
• Of course, it could be used. However, I doubt that it would be much fast, as you iterate only a few dozen times. Jan 11, 2018 at 18:52

FoldList won't give a big raise to the efficiency (the most time consuming part is equation solving), but it does make the code conciser:

sol[eta_] :=
ParametricNDSolveValue[{x''[t] + 1/10 x'[t] + x[t] + 1/100 x[t]^3 == 8/10 Cos[eta t],
x[0] == x0, x'[0] == v0,
WhenEvent[x'[t] < 0, Sow[x[t], eta], "DetectionMethod" -> "Interpolation"]},
{x[100], x'[100]}, {t, 0, 100}, {x0, v0}];
etalst = Range[8/10, 1, 1/100];
{iclst, maxlst} = Reap@FoldList[sol@#2 @@ # &, {0, 0}, etalst];

maxValue = Max /@ ({Partition[iclst\[Transpose] // First, 2, 1], maxlst}\[Transpose])

(* {2.77457, 3.50534, 2.82074, 2.40488, 2.96662, 3.81174, 4.24709, 4.08059, 3.62443,
2.8937, 3.46553, 4.14776, 4.67875, 4.97169, 4.68464, 3.92852, 4.15529, 4.71447, 5.13502,
5.53676, 5.53677} *)

• I figure my question was confusing, so I added some more code. Jan 12, 2018 at 10:51
• @freddy90 So you're mostly interested in speed, right? Then I think my method i.e. detecting inside NDSolve should be faster than yours, have you tried it? BTW, some other method for extracting maximum can be found here: mathematica.stackexchange.com/a/58256/1871 Jan 12, 2018 at 12:26
• Tried it out and you're right, it is faster than my approach and I learned a lot from your code. Thanks a lot!! Jan 16, 2018 at 10:02