Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]]
share|improve this question

2 Answers 2

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

share|improve this answer
    
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

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}]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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