1
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ My advice: Better use Do instead of For. Write it down as it comes to your mind and optimize later. $\endgroup$ Jan 11, 2018 at 17:54
  • $\begingroup$ 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? $\endgroup$
    – freddy90
    Jan 11, 2018 at 18:32
  • $\begingroup$ Of course, it could be used. However, I doubt that it would be much fast, as you iterate only a few dozen times. $\endgroup$ Jan 11, 2018 at 18:52

1 Answer 1

4
$\begingroup$

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} *)
$\endgroup$
3
  • $\begingroup$ I figure my question was confusing, so I added some more code. $\endgroup$
    – freddy90
    Jan 12, 2018 at 10:51
  • 1
    $\begingroup$ @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 $\endgroup$
    – xzczd
    Jan 12, 2018 at 12:26
  • $\begingroup$ Tried it out and you're right, it is faster than my approach and I learned a lot from your code. Thanks a lot!! $\endgroup$
    – freddy90
    Jan 16, 2018 at 10:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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