MATH529L Lab07: FFT data analysis

Forced linear oscillator

Consider point mass  affixed to one end of a spring of stiffness  with the other spring end fixed. Displacement from the equilibrium position due to external force f sin(ω t) is  and the equation of motion is

m y''+ k y = f sin(ω t)

Choose units of time, length, mass such that the equation is expressed in non-dimensional terms as

y''+  y =  sin(ω t)

DT=6π;
sol=TrigReduce[DSolve[{ y''[t]+y[t]==Sin[ω t],y[0]==0,y'[0]==0},y[t],t][[1,1]]]

The non-dimensional system eigenmodes of the unforced system have pulsation Ω=1=(Sqrt[k/m]), and the forced system exhibits resonance when . The solution is the superposition of the eigenmode frequency and the forcing frequency. No other frequency is present.

Plot[y[t] /. sol  /.  ω->2,{t,0,DT}]

Sampling the solution constructs a vector u

u=Table[ y[t] /. sol  /.  ω->2, {t,0,DT,π/1024}];

ListPlot[u,Joined->True]

Integrate[y[t] /. sol,{t,0,4π}]

m=Length[u]; m2=Floor[m/2]; U=Fourier[u]; P=Sqrt[U*Conjugate[U]]; ListLogPlot[P[[1;;m2]],PlotRange->All,Frame->True,GridLines->Automatic,FrameLabel->{"Pulsation number","Power"}]

The above power spectrum exhibits two peaks corresponding to the  pulsations.

Forced nonlinear oscillator

Consider now a point mass  affixed to one end of a spring of nonlinear stiffness ϵ  (y^3) with the other spring end fixed. Displacement from the equilibrium position due to external force f sin(ω t) is  and the equation of motion is

m y''+  k y + ϵ (y^3) = f sin(ω t)

Choose units of time, length, mass such that the equation is expressed in non-dimensional terms as

y''+  y + ϵ (y^3) =  sin(ω t)

solN=NDSolve[{ y''[t]+y[t]+ ϵ (((y[t]))^3)==Sin[ω t],y[0]==0,y'[0]==0} /. {ϵ->4×(10^(-2)),ω->2},y[t],{t,0,DT}][[1,1]]

The non-dimensional system eigenmodes of the unforced system have pulsation Ω=1=(Sqrt[k/m]), and the forced system exhibits resonance when . The solution is the superposition of the eigenmode frequency and the forcing frequency. No other frequency is present.

Plot[{y[t] /. sol  /.  ω->2, y[t] /. solN},{t,0,DT},PlotRange->All,Frame->True,GridLines->Automatic,FrameLabel->{"t","y"},PlotLabel->"Comparison of linear/nonlinear oscillators"]

Sampling the solution constructs a vector u

uN=Table[ y[t] /. solN , {t,0,DT,π/1024}];

ListPlot[{u,uN},Joined->True]

m16=Floor[m/16];UN=Fourier[uN]; PN=Sqrt[UN*Conjugate[UN]]; ListLogPlot[{P[[1;;m16]],PN[[1;;m16]]},PlotRange->{{1,24},All},Frame->True,GridLines->Automatic,FrameLabel->{"Pulsation number","Power spectra comparison"}]

The above power spectrum exhibits two peaks corresponding to the  pulsations.

Perturbation analysis

Consider  solution to

y''+  y + ϵ (y^3) =  sin(ω t)

as a power series in ϵ

y[t]= (Subscript[y, 0])[t]+ϵ (Subscript[y, 1])[t]+(ϵ^2) (Subscript[y, 2])[t]+..

DEq= y''[t]+y[t]+ ϵ (((y[t]))^3)==Sin[ω t]

pSer[n_]:= Sum[y[k,t] (ϵ^k),{k,0,n}]

pSer[2]

D[pSer[2],{t,2}]

nOrder=2;
DEqnOrder=DEq /. {y[t]->pSer[nOrder], y''[t]->D[pSer[nOrder],{t,2}]}

DEqnOrder /. ϵ->0

sol0=DSolve[{DEqnOrder /. ϵ->0,y[0,0]==0,Derivative[0,1][y][0,0]==0},y[0,t],t][[1,1]]

Coefficient[DEqnOrder[[1]],ϵ,1]

Coefficient[DEqnOrder[[2]],ϵ,1]

DEqHierarchy=Table[Coefficient[DEqnOrder[[1]],ϵ,k]==Coefficient[DEqnOrder[[2]],ϵ,k],{k,0,nOrder}]

sol0=TrigReduce[DSolve[{DEqHierarchy[[1]],y[0,0]==0,Derivative[0,1][y][0,0]==0},y[0,t],t][[1,1]]]

sol1=TrigReduce[DSolve[{DEqHierarchy[[2]] /. sol0,y[1,0]==0,Derivative[0,1][y][1,0]==0},y[1,t],t][[1,1]]]

sol2=TrigReduce[DSolve[{DEqHierarchy[[3]] /. sol0 /. sol1,y[2,0]==0,Derivative[0,1][y][2,0]==0},y[2,t],t][[1,1]]]

Multidimensional, experimental power spectra

Experiments often furnish data on systems, such as the video recordings of capillary gravity waves under external mono-chromatic excitation.

vid1=Import["/Users/mitran/courses/MATH529L/GravCapWave1.mp4"]

im=VideoExtractFrames[vid1,76]

ImageQ[im]

ImageDimensions[im]

imBW=ColorConvert[im,"Grayscale"]

window=ImageTake[imBW,{501,1200},{141,810}]

A=ImageData[window]; Dimensions[A]

ListPlot[A[[200,90;;580]]]

xSig=A[[200,90;;580]];

Xsig=Fourier[xSig]; PXsig=Sqrt[Xsig*Conjugate[Xsig]]; m=Length[PXsig]; m4=Floor[m/4];
ListLogPlot[PXsig[[1;;m4]],Frame->True,GridLines->Automatic,FrameLabel->{"kx","Power"}]

vid2=Import["/Users/mitran/courses/MATH529L/GravCapWave2.mp4"]

Length[vid1]

Tasks

1. Compare numerical solution of nonlinear oscillator to perturbation solution.

2. Construct power spectra along  of gravity-capillary waves in the close to linear regime.

3. Construct power spectra along  of gravity-capillary waves in the non-linear regime.