A question elsewhere on this site asks about detecting the frequency of oscillations in a system defined by differential equations.

The equation is $y'(x) = y(x) \cdot \cos(x + y(x))$. The solution can be found numerically (plot of $y(x)$ over $x$):

enter image description here

In my answer, I reasoned that the period may be found by taking one of the peaks on an FFT. However, the following Mathematica code for taking a Fourier transform of a bunch of samples seems to give nonsense results:

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 10000 equally spaced points *)
numPoints = 10000;
samplePoints = Range[0, 30, 30/numPoints];
samples = Evaluate[y[samplePoints] /. s][[1]];

ft = Fourier[samples];
ListPlot[{Abs[ft], Re[ft], Im[ft]}]

Result (blue is absolute, purple is real part, yellow is imaginary part):

enter image description here

Since the x-axis is supposed to represent frequency, I can also convert it to period:

t = 1/Fourier[samples];
ListPlot[{Abs[t], Re[t], Im[t]}]

Which still fails to produce a noticeable signal which would correspond to the real period of 6.44:

enter image description here

The problem is probably not related to Mathematica (because Matlab produces similar results), but has to do with my poor understanding of the FT.

Conceptually, I think it is reasonable to expect the FT to produce a strong signal at the frequency corresponding to the real period of the plot, and its harmonics.

Why then, is this input producing a "spectrum" composed of complex numbers (specifically, what is the intuitive/practical interpretation of a complex result)? Why am I unable to interpret the result of the FT to get the correct answer?


1 Answer 1


This solution is not periodic. Instead, the oscillations are decaying. In fact, to second order in $y$ we have $$ y' = \cos(x) y - \sin(x) y^2 $$ which has solutions $$ y = \dfrac{e^{\sin(x)}}{C + \int e^{\sin(x)}\sin(x)\; dx}$$ Now $\int_{0}^{2\pi} e^{\sin(x)} \sin(x)\; dx > 0$, call this $2 \pi r$ (in fact $r = I_1(1)$ where $I_1$ is a modified Bessel function). Then $$\int_0^x e^{\sin(x)} \sin(x)\; dx = r x + O(1)$$ so that (for this second order approximation) $$ y = \dfrac{e^{\sin(x)}}{r x} + O(1/x^2)$$ I would expect that for the original ODE we have something similar.
In fact here is a plot of $x y$ as a function of $x$ for a numerical solution with initial value $y(0) = 1$:

enter image description here

So if you're looking for periodicity, you might look not at $y$ but at $xy$ (after waiting some time for the solution to settle down).

  • $\begingroup$ Nice, and the approximation is a good match to the original DE too. $\endgroup$ May 30, 2014 at 21:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .