0

I have a graph here and I want to be able to manipulate the 4 parameters that are in the Block command. I want to watch the graph change as I change the parameters.

I've done this before, but the Block command that Mathematica suggested throws me off a little bit.

Here is the code:

 test = Block[{ε = 3, ω0 = 1, ωf = 1.788, 
 A = 5}, Reap[
 NDSolve[{x''[t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] + 
     A*Sin[ωf*t], x[0] == 1, x'[0] == 0, 
   WhenEvent[Mod[t, 2 π/ωf] == 0, 
    Sow[{x[t], x'[t]}]]}, {}, {t, 0, 10000}, 
  MaxSteps -> ∞]]][[-1, 1]];

ListPlot[test, ImageSize -> Medium, PlotRange -> All, 
 PlotStyle -> PointSize[0.0035]]
Corey Kelly
  • 1,738
  • 9
  • 23
Slightly
  • 321
  • 1
  • 2
  • 10

2 Answers2

1

Memoize NDSolve

If you want to convince with a Manipulation it must be smooth. But I'm sorry to say, that if you naively choose all the parameters and use them for Manipulate this is frustrating slow.

The thing I'm trying here is to memoize the output of NDSolve inside a DynamicModule, which lives inside a Manipulate.

This is the first time I've done this and it is actually working pretty well (with lots of help of Mma I guess...).

Manipulate[
    DynamicModule[{},
        ListPlot[lpFunc[eps, om0, omf, A][[-1, 1]], ImageSize -> Medium, 
                 PlotRange -> All, PlotStyle -> PointSize[0.0035]],
        Initialization :> (lpFunc[eps_, om0_, omf_, a_] := 
           lpFunc[eps, om0, omf, a] = 
              Reap[NDSolve[{x''[t] == 
                eps (1 - x[t]^2) x'[t] - om0^2*x[t] + A*Sin[omf*t], 
                x[0] == 1, x'[0] == 0, 
                WhenEvent[Mod[t, 2 π/omf] == 0, 
                Sow[{x[t], x'[t]}]]}, {}, {t, 0, 10000}, 
                MaxSteps -> ∞]]])
   ],
   {{eps, 3, ε}, 1, 5},
   {{om0, 1, ω0}, 1, 5},
   {{omf, 1.788, ωf}, 0, 5},
   {{A, 5}, 1, 10},
   TrackedSymbols -> {eps}
]

Your differential equation has very easily some singularities, so you better limit the range of what a user may choose, but this is your choice. I've chosen eps as the trigger to re-evaluate the plot (if it isn't already memoized).

Update: Gain some speed!

Although I think that the solution so far is nice, it is still not a decent Manipulate experience. So not a good advertisement.

What I'd do if I were you is the following:

1) Rewrite the function for NDSolve: (please have a look first at the bottom of this post, where @rcollyer rewrote the function inside a Compile block, gaining a massive speed improvement.)

lpFunc[eps_, om0_, omf_, A_] := 
lpFunc[eps, om0, omf, A] = 
    Reap[NDSolve[{x''[t] == eps (1 - x[t]^2) x'[t] - om0^2*x[t] + A*Sin[omf*t], x[0] == 1, 
        x'[0] == 0, WhenEvent[Mod[t, 2 Pi/omf] == 0, Sow[{x[t], x'[t]}]]}, {}, {t, 10000},
        MaxSteps -> Infinity, 
        Method -> {Automatic, "SymbolicProcessing" -> 0}
    ]
]

I wanted to make NDSolve faster, that's why I turned off SymbolicProcessing, since this will be evaluated by the optimizer hundred of times.

2) Generate a table of precalculated results:

ParallelTable[Quiet@lpFunc[eps, om0, omf, a], {eps, 1, 5, 0.2}, {om0, 1, 5, 0.2}, 
    {omf, 0, 5, 0.2}, {a, 1, 10, 0.2}];

Now it's time to go out for lunch break or anything else, but it must take time...

3) DumpSave the symbol lpFunc to a .mx file:

SetDirectory[$TemporaryDirectory];
DumpSave["preprocessed.mx", lpFunc];

4) Create a new notebook and import the .mx file and define your Manipulate:

SetDirectory[$TemporaryDirectory];
<< preprocessed.mx

5) Create your Manipulate:

Manipulate[
    ListPlot[lpFunc[eps, om0, omf, A][[-1, 1]], ImageSize -> Medium, 
        PlotRange -> All, PlotStyle -> PointSize[0.0035]],
    {{eps, 3., ε}, 1, 5},
    {{om0, 1., ω0}, 1, 5},
    {{omf, 1.7, ω0f}, 0, 5},
    {{A, 5.}, 1, 10}
]

This will be now a nice Manipulate experience, where your plot will change smoothly. Something you can show.

Update End

By the way. I've had ancient greek at school and it gives me bad memories if someone is using greek letters (although not as bad as the memories for latin...phew)

Please have a look at this post where some of the authorities here made decent effort on importing unicode letters. Another nice hint from @cormullion...What I'd do without him...

I hope this helps or at least gives an idea for what you try to achieve.

Update 2: Gain Extra Speed in Prep

Since the function takes numbers and produces numbers, you can compile it as a numeric function into C. Which makes it FAR faster. Do it so:

cFunc = 
  Compile[{eps, om0, omf, A}, 
   Reap[NDSolve[{x''[t] == 
       eps (1 - x[t]^2) x'[t] - om0^2*x[t] + A*Sin[omf*t], x[0] == 1, 
      x'[0] == 0, 
      WhenEvent[Mod[t, 2 Pi/omf] == 0, Sow[{x[t], x'[t]}]]}, {}, {t, 
      10000}, MaxSteps -> Infinity, 
     Method -> {Automatic, "SymbolicProcessing" -> 0}]], 
   CompilationTarget -> "C"];
lpFunc[eps_, om0_, omf_, A_] := lpFunc[eps, om0, omf, A] = cFunc[eps, om0, omf, A];

I'm seeing a 2-3X improvement in speed of original data generation.

End Update

Stefan
  • 5,347
  • 25
  • 32
  • Re: greek; also see this – cormullion Jun 26 '13 at 16:58
  • Thanks :) didn't knew they've made progress on that matter. I failed utterly porting @Mr.Wizard's entry to use pbcopy (for Mac), so I gave up on Unicode. Will check it, although there are bad memories concerning Homer and Hexameter ;) – Stefan Jun 26 '13 at 17:14
  • @rcollyer your edit is simply amazing. thank you for that...totally lost Compile out of sight when i wrote this. thanks a lot. now this post seems to be pretty decent. what do you think? – Stefan Jun 26 '13 at 20:50
  • Hey, if it was my doing, I'd gladly take credit. But, alas and alack, I merely fixed up someone else's stuff. Blame them. :) – rcollyer Jun 26 '13 at 21:05
  • yep, sorry, guilty as charged. :-) If I see a good place to plug in Parallelize or Compile, I can't help myself. :-) – Gregory Klopper Jun 26 '13 at 23:38
  • @GregoryKlopper NDSolve is not compilable. I didn't try it, but are you really sure you see a performance improvement, given that the compiled function immediately calls back out of the VM to evaluate the top-level code? If you've been experimenting, you might easily have some memoized definitions still hanging around, and thus fool yourself over its actual performance. – Oleksandr R. Jun 27 '13 at 01:49
  • When I memorize the definitions, the answer gets back in <<1sec. When I do not, the difference in time is 2-3 fold, but still 10-30 seconds. Not sure how to explain it. Wonder if CompilationTarget->"C" has anything to do with it? – Gregory Klopper Jun 27 '13 at 03:49
  • @GregoryKlopper ok. so if you found out compile does really help here, or it was just the memoized effect we might update this post.... – Stefan Jun 27 '13 at 06:01
  • My tests show, as I commented above, that it really does have an effect on speed. As any good science experiment, it should be repeatable by others. I invite someone to copy-paste the code, and confirm their timings. I used Remove@lpFunc before every re-assignment of lpFunc to do the testing, but I don't remember if I killed the kernel or didn't. Shouldn't matter. One can hardly mistaken 15 seconds for 0.01s, can they? :-) – Gregory Klopper Jun 27 '13 at 22:51
0

Hope this is what you wanted to accomplish. Good luck!

Manipulate[
 test = Reap[
    NDSolve[{x''[
        t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] + 
        A*Sin[ωf*t], x[0] == 1, x'[0] == 0, 
      WhenEvent[Mod[t, 2 π/ωf] == 0, 
       Sow[{x[t], x'[t]}]]}, {}, {t, 0, 10000}, 
     MaxSteps -> ∞]][[-1, 1]];
 ListPlot[test, ImageSize -> Medium, PlotRange -> All, 
  PlotStyle -> PointSize[0.0035]]
 , {{ε, 3}, 1, 5}
 , {{ω0, 1}, 1, 5}
 , {{ωf, 1.788}, 0, 5}
 , {{A, 5}, 1, 10}]
Corey Kelly
  • 1,738
  • 9
  • 23
Gregory Klopper
  • 1,370
  • 9
  • 21