# Differential Equation power spectra

ε = .1; A = 0; ω0 = 1; ωf = 0;
data =
NDSolve[
{x''[t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] + A*Cos[ωf*t], x[0] == 1, x'[0] == 0},
{t, 0, 100}];
Periodogram[data, ScalingFunctions -> "Absolute", PlotRange -> All]


How would I create a suitable Power spectrum of the frequencies that this system has. For these set of values (at the top), there should ideally be only one unique frequency (ω0) since the other one is not in play.

Mathematica's example of power spectra utilizes the periodogram funcition, although I don't really understand what it does mathematically, so there is room for using another type of graph.

ε = 1/10; A = 0; ω0 = 1; ωf = 0;
sol = NDSolve[{x''[
t] == ε*(1 - x[t]^2) x'[t] - ω0^2*x[t] +
A*Cos[ωf*t], x[0] == 1, x'[0] == 0},
x[t], {t, -100, 100}];
Plot[x[t] /. sol, {t, -100, 100}]


Periodogram[Flatten[Table[x[t] /. sol, {t, -100, 100, 0.5}]],PlotRange -> All]


Here's my shot at it:

ε = .1; A = 0; ω0 = 1; ωf = 0;

sol = First[NDSolve[
{x''[t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] + A*Cos[ωf*t], x[0] == 1, x'[0] == 0},
x, {t, 0, 100}]];
f = x /. sol;
data = Table[f[t], {t, 0, 100}];

Periodogram[data, ScalingFunctions -> "Absolute"]


• I've worked with this one, but it does not give me the correct results. – Slightly Jul 3 '13 at 13:54
• A-ha... And we should guess: Are there any error messages ? What do they say ? Any attempts to fix the "problems" ? – Sektor Jul 3 '13 at 14:06
• There are no error messages, the problem is that I am not getting a frequency on the x-axis, and the situations where I expect chaos, the power spectrum should be more messy. – Slightly Jul 3 '13 at 19:56
• Are you using the same code ? If not, please, do provide the code you are using. – Sektor Jul 4 '13 at 17:07
• Well, I've been trying it with yours, but there is no luck, and I tried with the other answer, but everything should be above the x axis. – Slightly Jul 4 '13 at 21:33

Consider the equation: $$\ddot{x} = (x,\dot{x},t)$$ The numerical solution generates a sequence of values of x at discrete values of t. This sequence is known as time series. Plotting x pk t might seem random or noisy, but often it contains harmonics at certain frequencies. Dominant frequencies in time series data can be investigated using the Discrete Fourier transform. For example we take N equally spaced values from the time series - $$x = x_{0},x_{1}, ..., x_{N}$$ The Discrete Fourier transform is defined as follows: $$X_{k}=\frac{1}{\sqrt{N}}\sum_{m=0}^{N}x_{m}e^{\frac{-2\pi ikm}{N}}, (k=0,1,2,...,N-1)$$ Generally $X_{k}$ is a complex number. To analyse the frequency structure we analyse the power spectrum $P(\omega _{k})$ defined by: $$P(\omega _{k})=X_{k}\overline{X_{k}}=\left | X_{k} \right |^{2}$$

ε = .1; A = 0; ω0 = 1; ωf = 0;

ff = x[t] /.
First[NDSolve[{x''[
t] == ε (1 - x[t]^2) x'[t] - ω0^2*x[t] +
A*Cos[ωf*t], x[0] == 1, x'[0] == 0},
x, {t, 0, 300}]][[1]]

Table[ff, {t, 0, 100}] // Fourier // Abs //
ListLogPlot[#, PlotRange -> All] &


• I feel as if we are getting closer to the answer as I see you are using Fourier, although this isn't getting me the correct results. Here, there should only be one frequency at around .1. Can you tell me what you are plotting on the x and y axes? – Slightly Jul 19 '13 at 15:32
• @Slightly Correct me if I am wrong, but this is the answer to your question. Nothing more, nothing less. – Sektor Jul 22 '13 at 13:19
• This is getting closer to what I am looking for, although the values are not correct. Since I have ω0 = 1, I should only have one peak at 1 (on the x axis) – Slightly Jul 23 '13 at 2:10
• I don't agree. Thanks for the help, though. – Slightly Jul 23 '13 at 15:48

The code that gives me the best results is this.

ε = 3; A = 5; ω0 = 1; ωf = 1.788;
sol = NDSolve[{x''[t] == ε*(1 - x[t]^2) x'[t] - ω0^2*x[t] +
A*Sin[ωf*t], x[0] == 1, x'[0] == 0}, x[t], {t, 0, 500}, MaxSteps -> Infinity]

Periodogram[Flatten[Table[x[t] /. sol, {t, 0, 500, 1}]],
PlotRange -> All, ScalingFunctions -> "Absolute"]

• Well, you could have mentioned the parameters are different ... – Sektor Jul 23 '13 at 19:35
• It's a completely different graph. – Slightly Jul 23 '13 at 22:44
• Yeah, because the conditions are different :D :D + There is a Sin[$\omega ft$] term in your answer whereas there's a Cos[$\omega ft$] in our code, based on the equation you provided in the first place. Obviously you made a lot of mistakes and you did not acknowledge making them. – Sektor Jul 24 '13 at 7:40
• Oh yes, I do agree that the forcing function is different, but it does not change anything in the output. – Slightly Jul 24 '13 at 13:23