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 am wondering how I could determine the frequency of oscillations of a differential model equation? How could I find the frequency from this example given in Mathematica Documentation:

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}]
Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]

Mathematica graphics

I am asking because I need to determine the frequency of one model term as a function of a the strength of a second model term, in a much more complicated model than described above.

share|improve this question
up vote 17 down vote accepted

I'd advocate taking differences between successive peaks and likewise troughs. These can be found by keeping track of when the derivative is zero.

pts = 
 Reap[s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1, 
      WhenEvent[y'[x] == 0, Sow[x]]}, {y, y'}, {x, 0, 30}]][[2, 1]]

(* Out[290]= {0.448211158984, 4.6399193764, 7.44068279785, 10.953122261, \
13.8722260952, 17.2486864443, 20.2244048853, 23.5386505821, \
26.5478466115, 29.8261176372} *)

Plot[{Evaluate[y[x] /. s], Evaluate[y'[x] /. s]}, {x, 0, 30}, 
 PlotRange -> All]

enter image description here

diffs = Differences[pts, 1, 2]

(* Out[288]= {6.99247163887, 6.31320288463, 6.43154329733, \
6.29556418327, 6.35217879014, 6.28996413777, 6.32344172616, \
6.28746705515} *)


(* Out[289]= 6.41072921417 *)

Noticing that the first might be an outlier, one might choose to discard it.


(* Out[291]= 6.32762315349 *)
share|improve this answer
Thanks, that works perfectly! – tarhawk May 15 '14 at 19:47

The classical Fourier approach

Your data :

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}][[1, 1, 2]];
Plot[s[x], {x, 0, 30}, PlotRange -> All]

enter image description here

The horizontal axis is supposed to be graduated in seconds.
First, we sample s at the frequency 5 Hz :

data00 = s /@ Range[0, 30, 1/5];
ListLinePlot[data00, PlotRange -> All, PlotStyle -> Thick, Joined -> False]

enter image description here

To increase the resolution of the Fourier transform, this record is extended to 1000 samples, so that the resolution of the fourier transform is fine enough (5/1000 = 5 mHz between each sample) :

data01 = PadRight[data00, 1000];
ListLinePlot[data01, PlotRange -> All, PlotStyle -> Thick]

enter image description here

The spectrum is calculated :

data02 = Abs[Fourier[data01]];
ListLinePlot[data02, PlotRange -> All]

enter image description here

To localise precisely where is the peak near the coordinates {40, 0.45}, we reject the uninteresting frequencies. Then we look for the maximum with Ordering[] :

data03 = Table[If[20 < i < 500, data02[[i]], 0.], {i, Length[data02]}];
ListLinePlot[data03, PlotRange -> All] 
First[Ordering[data03, -1]] (* maximum *)

enter image description here


The peak is at position 32 on the horizontal axis. This means that the frequency is (32 - 1)*0.005 = 0.155Hz, that is to say period=6.45 second.

share|improve this answer
Beat me to it! Good use of Fourier. – blochwave May 16 '14 at 17:43

I'm not very good at Mathematica, but I think there are two "usual" ways of solving this problem.

One is to do a Fourier transform of your data, and find the peak. Unfortunately I couldn't get a result from Mathematica in a way that makes sense for this.

The other option is to use an autocorrelation plot, and take the best peak. I'm not sure how to do this for a function, but you can sample points along your function and then work with those points:

Let's input your data, and make a plot: (I won't include a plot here since you already have it in the question)

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}]; 
Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]

Sample the solution at 10.000 equally spaced points:

numPoints = 10000;
samplePoints = Range[0, 30, 30/numPoints];
samples = Evaluate[y[samplePoints] /. s][[1]];
ListPlot[samples]; (* You can plot the sampled points to be sure, but I disabled it *)

Calculate the autocorrelation:

autocorrelations = ListCorrelate[samples, samples , {1, 1}, 0];

Find peaks by looking for points where the differential changes sign:

peaks = 1 + Position[Differences[Sign[Differences[autocorrelations]]], -2];
bestPeriods = Extract[samplePoints, peaks];

Note that if you have Mathematica 9+, then the above code can be simplified to just MinDetect. See

Create a nice plot of the autocorrelation, and draw a red line at the most likely peak:

xyPoints = Transpose[{samplePoints, autocorrelations}];
ListLinePlot[xyPoints, PlotRange -> All, GridLines -> {{{First[bestPeriods], Red}}, {}}, PlotLabel -> "Autocorrelation plot - the most likely period is " <> ToString[First[bestPeriods]] <> ":"]

enter image description here

If you want to use the detected period for something, you can access it with First[bestPeriods].

share|improve this answer
This is nice. I also would like to see an approach using Fourier. I think it can be done but was also not able to get something that worked. – Daniel Lichtblau May 16 '14 at 14:42

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.