Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have numerically solved and plotted the solution to the Mathieu equation using NDSolve successfully, as below:

a = 0.5;
b = 0.1;
w0 = 2.0;
T = 2 Pi/w0;

f[t_] := f[t] = y''[t] + (a + b  Cos[t]) y[t] == 0

s = NDSolve[{f[t], y[0] == 1, y'[0] == 0}, y, {t, 0, 1}]

Plot[Evaluate[y[t] /. s], {t, 0, 1}, PlotRange -> All]

But now I wish the normalize the period so I can plot from t=0 -> t=1 and get 1 full period. I can do this by hand, but how do I do this inside of the NDSolve function?

share|improve this question
What are w0 and T? You don't seem to be using them ... – Dr. belisarius Mar 25 '14 at 17:42
Mathematica knows MathieuC. – Kuba Mar 25 '14 at 17:47
w0 is the normal frequency of the periodic function, in this case Cosine. T would be the normalization factor, but if I multiply the Cos function by T, my plot doesn't normalize as I expected. In fact, nothing changes. – gKirkland Mar 25 '14 at 17:49
@Kuba Yes, I know, but I will use this code to check more complex DE's with periodic coefficients later so I'd like to make it general instead of using MathieuC. – gKirkland Mar 25 '14 at 17:51
I see a difference with and without the factor T. Look at the scale on the vertical axis. – Michael E2 Mar 25 '14 at 18:27
up vote 4 down vote accepted

Is this what you are after? Using WhenEvent to find the extrema..

 s = Reap[ 
      NDSolve[{y''[t] + (a + b Cos[t]) y[t] == 0, y[0] == 1, y'[0] == 0, 
         WhenEvent[ y'[t] == 0 && y[t] > 0 , Sow[{ t, y[t]}]]}, 
        y, {t, 0, 50}] ]
   Plot[Evaluate[y[t] /. s[[1]]], {t, 0, 50 }, PlotRange -> All],
      ListPlot[s[[2, 1]], PlotStyle -> {Red, PointSize[.02]}]

enter image description here

The solution isn't exactly periodic, or it has a long its not exactly clear what you mean to normalize to the period.

If I put

WhenEvent[ y'[t] == 0 && Abs[y[t] - 1] < .001 , Sow[{ t, y[t]}]]

we find "period" of ~ 170 where the solution returns to a good approximaiton to y=1,y'=0

share|improve this answer

Your Answer


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.