# Plot FFT of an interpolating function

Posted 9 years ago
7834 Views
|
5 Replies
|
9 Total Likes
|
 Hi everyone, I have an oscillating system and i have to find the spectral components of motion of each oscillator. I found xn(t) solving system with NDSolve, I tried to calculate FFT on that result, but i cannot plot the FFT result. Could you help me please to understand what I'm doing wrong? Thank you in advance, here's the code:  sistema = {m*x1''[t] == 0, m*x2''[t] == 0, m*x3''[t] == 0, m*x4''[t] + k1*(-x2[t] - x6[t] + 2 x4[t]) == 0, m*x6''[t] + k1*(-x5[t] - x4[t] + 2 x6[t]) == 0, m*x5''[t] + k1*(-x3[t] - x6[t] + 2 x5[t]) == 0, m*y1''[t] == 0, m*y2''[t] == 0, m*y3''[t] == 0, m*y4''[t] + k1*(-y2[t] - y6[t] + 2 y4[t]) == 0, m*y6''[t] + k1*(-y5[t] - y4[t] + 2 y6[t]) == 0, m*y5''[t] + k1*(-y3[t] - y6[t] + 2 y5[t]) == 0, m*z1''[t] == 0, m*z2''[t] == 0, m*z3''[t] == 0, m*z4''[t] + k1*(-z2[t] - z6[t] + 2 z4[t]) == 0, m*z6''[t] + k1*(-z5[t] - z4[t] + 2 z6[t]) == 0, m*z5''[t] + k1*(-z3[t] - z6[t] + 2 z5[t]) == 0, x1'[0] == 0, x2'[0] == 0, x3'[0] == 0, x4'[0] == 0, x5'[0] == 0, x6'[0] == 0, y1'[0] == 0, y2'[0] == 0, y3'[0] == 0, y4'[0] == 0, y5'[0] == 0, y6'[0] == 0, z1'[0] == 0, z2'[0] == 0, z3'[0] == 0, z4'[0] == 0, z5'[0] == 0, z6'[0] == 0, x1[0] == 5.62726, x2[0] == 5.71947, x3[0] == 6.78594, x4[0] == 6.96736, x5[0] == 8.03382, x6[0] == 8.12604, y1[0] == 3.23074, y2[0] == 4.52119, y3[0] == 2.41941, y4[0] == 4.99979, y5[0] == 2.89802, y6[0] == 4.18846, z1[0] == 1.12774, z2[0] == 0.5448, z3[0] == 1.23874, z4[0] == 0.07373, z5[0] == 0.76767, z6[0] == 0.18473}; soluzione = NDSolve[sistema, {x1[t], x2[t], x3[t], x4[t], x5[t], x6[t], y1[t], y2[t], y3[t], y4[t], y5[t], y6[t], z1[t], z2[t], z3[t], z4[t], z5[t], z6[t]}, {t, 0, 1}, MaxSteps -> 1000000] fun = x4[t] /. soluzione; pointNo = 2560000000000; inc = 1/pointNo; data = Table[fun, {t, 0., 2.7*10^-10, inc}]; spettro = Fourier[data] ListLinePlot[Abs[Fourier[data]], PlotRange -> All] The figure I obtain is a stack of points.
5 Replies
Sort By:
Posted 9 years ago
 Forget to mention that NDSolveValue[] can be used to generate data instead of WhenEvent[]. The latter is probably more suitable for more sophisticated ''events''. The modified code is a little bit cleaner and tiny bit faster: DATA = NDSolveValue[ { x1'[s] == y1[s], y1'[s] == - 16. x1[s], x2'[s] == y2[s], y2'[s] == - 25. x2[s], x1[0] == 0.0, y1[0] == 1.0, x2[0] == 1.0, y2[0] == 0.0 }, {x1[#],x2[#]} & /@ Table[s,{s,1.,1000.,1.}], {s,0,1000}, MaxSteps -> Infinity, Method -> { "ImplicitRungeKutta", "DifferenceOrder" -> 10, "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", "ImplicitSolver" -> { "Newton", AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, "IterationSafetyFactor" -> 1 } } ] // Transpose ; ListLogPlot[Abs[Fourier/@DATA],PlotRange->All,Frame->True,ImageSize-> 500] I.M.
Posted 9 years ago
 Consider this new feature in Mathematica if you want to include {m, k1} as a set of variables for your ODE's.
Posted 9 years ago
 Hi, Phil,Fourier spectra are supposed to be different for different sample length. Read about fourier frequency resolution and Nyquist frequency to convert x axis into frequency. Also, take a look at Fourier[] documentation, FourierParameters and Applications section in particular.I.M.
Posted 9 years ago
 Hi Ivan, k1 and m are two numeric variables. I tried to modify my code as you suggested and now it seems to work. That example was very useful. Thank you so much. There's only an issue: why if I change s step, the result is totally different? For example i tried to change in your code [Mod[s,0.1] == 0. and the new plot that I obtain is that: which is totally different from that: How values on the x axis are related to the frequencies of the system?
Posted 9 years ago
 Hi, Phil,Not sure what is wrong with your code (you forget to include definitions for m and k1). Here is a simple example for 2 uncoupled harmonic oscillators: DATA = Reap[ NDSolve[ { x1'[s] == y1[s], y1'[s] == - 16. x1[s], x2'[s] == y2[s], y2'[s] == - 25. x2[s], x1[0] == 0.0, y1[0] == 1.0, x2[0] == 1.0, y2[0] == 0.0, WhenEvent[Mod[s,1.] == 0.,Sow[{x1[s],x2[s]}]] }, {}, {s,0,1000}, MaxSteps -> Infinity, Method -> { "ImplicitRungeKutta", "DifferenceOrder" -> 10, "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", "ImplicitSolver" -> { "Newton", AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, "IterationSafetyFactor" -> 1 } } ] ][[2,1]] // Transpose ; ListLogPlot[Abs[Fourier/@DATA],PlotRange->All,Frame->True,ImageSize-> 500] Note, that I do not use interpolating function, but WhenEvent[] function to generate trajectory sample at s = 1,2, ... I.M.