Message Boards

WOLFRAM COMMUNITY

9623 Views

5 Replies

9 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Plot FFT of an interpolating function

Phil R

Posted 10 years ago

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 = {mx1''[t] == 0, mx2''[t] == 0, mx3''[t] == 0, mx4''[t] + k1(-x2[t] - x6[t] + 2 x4[t]) == 0, mx6''[t] + k1(-x5[t] - x4[t] + 2 x6[t]) == 0, mx5''[t] + k1(-x3[t] - x6[t] + 2 x5[t]) == 0, my1''[t] == 0, my2''[t] == 0, my3''[t] == 0, my4''[t] + k1(-y2[t] - y6[t] + 2 y4[t]) == 0, my6''[t] + k1(-y5[t] - y4[t] + 2 y6[t]) == 0, my5''[t] + k1(-y3[t] - y6[t] + 2 y5[t]) == 0, mz1''[t] == 0, mz2''[t] == 0, mz3''[t] == 0, mz4''[t] + k1(-z2[t] - z6[t] + 2 z4[t]) == 0, mz6''[t] + k1(-z5[t] - z4[t] + 2 z6[t]) == 0, mz5''[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.710^-10, inc}]; spettro = Fourier[data] ListLinePlot[Abs[Fourier[data]], PlotRange -> All] The figure I obtain is a stack of points.

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. enter image description here

POSTED BY: Phil R

5 Replies

Sort By:

Ivan Morozov

Ivan Morozov, Budker Institute of Nuclear Physics

Posted 10 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.

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.

POSTED BY: Ivan Morozov

Phil R

Posted 10 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 BY: Phil R

Ivan Morozov

Ivan Morozov, Budker Institute of Nuclear Physics

Posted 10 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 BY: Ivan Morozov

Shenghui Yang

Shenghui Yang, WOLFRAM

Posted 10 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 BY: Shenghui Yang

Ivan Morozov

Ivan Morozov, Budker Institute of Nuclear Physics

Posted 10 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.

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 BY: Ivan Morozov

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback