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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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