Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
ε = .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.

share|improve this question
add comment

4 Answers

ε = 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}]

function plot

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

periodogram plot

share|improve this answer
add comment

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"]

enter image description here

share|improve this answer
    
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
show 8 more comments

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] &

Plot

share|improve this answer
    
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 think you understand what I wrote above and what you actually need. To analyse the spectrum you need to construct $P( \omega_{k})$. As such it takes $\omega_{k}$ as a parameter and then produces the respective $P( \omega_{k})$. In your example $\omega_{0}$ is nothing more than a parameter you pass to the function. + "the only one peak" is generally speaking incorrect, because the power spectrum plot is symmetrical. This completely answers your question, you can accept an answer. –  Sektor Jul 23 '13 at 8:45
    
I don't agree. Thanks for the help, though. –  Slightly Jul 23 '13 at 15:48
add comment
up vote -3 down vote accepted

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"]
share|improve this answer
1  
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
1  
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
add comment

Your Answer

 
discard

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.