Resonance Phenomena

Resonance phenomena simulated by finite differences and biased noise inputs

In Lecture-21.nb, a second-order differencing scheme that iteratively solved  y'' + βy' + γy=0  with the specification of two initial values.
Modify this to add a little random noise to y[i] at each step and observe how this behaves---this version will store the noise added at each iteration so that it can be visualized later....

GrowListGeneralNoise[ValuesList_List,Δ_,α_,β_,randomamp_]:=Module[{Minus1=ValuesList[[1,-1]],Minus2=ValuesList[[1,-2]],noise=RandomReal[{-randomamp,randomamp}]},{Append[ValuesList[[1]],2 Minus1-Minus2+Δ (β (Minus2-Minus1)-α Δ Minus2)+noise],Append[ValuesList[[2]],noise]}]

Setting up a function that takes particular parameters and a "noise amplitude" of (10^(-5))

GrowListSpecificNoise[InitialList_List]  := GrowListGeneralNoise[InitialList,.001,2,0,10^(-5)]

Nest[GrowListSpecificNoise,{{1,1},{0,0}},10]

TheData = Nest[GrowListSpecificNoise,{{1,1},{0,0}},20000];

ListPlot[TheData[[1]]]

ListPlot[TheData[[2]]]

Now suppose there is a periodic bias  that tends to kick the displacement one direction more than the other:

GrowListBiasedNoise[ValuesList_List,Δ_,α_,β_,randomamp_,lambda_]:=Module[{Minus1=ValuesList[[1,-1]],Minus2=ValuesList[[1,-2]],biasednoise=0.5` randomamp (Cos[((2 π Length[ValuesList[[1]]])/lambda)]+RandomReal[{-1,1}])},{Append[ValuesList[[1]],2 Minus1-Minus2+Δ (β (Minus2-Minus1)-α Δ Minus2)+biasednoise],Append[ValuesList[[2]],biasednoise]}]

GrowListSpecificBiasedNoise[InitialList_List]  := GrowListBiasedNoise[InitialList,.001,2,0,10^(-6),4500]

Generate the data set---this takes quite a while

TheBiasedData = Nest[GrowListSpecificBiasedNoise,{{1,1},{0,0}},20000];

ListPlot[TheBiasedData[[1]]]
ListPlot[TheBiasedData[[2]]]

Resonance Phenomena by Solution of the ODE

The "periodic forcing function" can be any periodic function (i.e., a fourier series), but it is a bit simpler to analyze the effect of each mode separately.  Below, the forcing function will be assumed to be (Subscript[F, app])cos((Subscript[ω, app])t)
Solve problems in terms of the mass and natural frequency--eliminate the spring constant in equations by defining it in terms of the mass and natural frequency.

Kspring =  M ((ωchar)^2)

Mathematica can solve the nonhomogeneous ODE with a  forcing function at with an applied frequency:

yGeneralSol = Simplify[y[t]/.DSolve[M y''[t] + η y'[t] + Kspring y[t] == Fapp Cos[ωapp t],y[t],t][[1]]]

Consider the form of the above solution: The homogeneous part of the solution  (i..e, the part that depends on constants C[1] and C[2]) depends only on  the characteristic frequency, while the heterogeneous part includes the forcing and the characteristic frequency.   From a first glance at the denominator in the heterogeneous term, it appears the amplitudes are  maximized as ωapp ⟶ ωchar; and , it appears that,  as long as the viscous term η is non-zero, the solution is bounded. However, this is not the case as will be shown below.

Consider the behavior of the general solution at time t=0.  This will show that the homogeneous parts of the solution are needed to satisfy boundary conditions, even if the oscillator is initially at rest at zero displacement (i.e., y(0) = (OverDot[y])(0) = 0).

Simplify[yGeneralSol/.t->0]

Therefore,  the C[1] and C[2] in the general solution, must pick up terms that depend on  ωapp and these must be included in considerations of the limit ωapp ⟶ ωchar.

Consider the particular case of an oscillator at rest at zero displacement:

yParticularSol = Simplify[y[t]/.DSolve[{M y''[t] + η y'[t] + Kspring y[t] == Fapp Cos[ωapp t], y[0]==0, y'[0]==0},y[t],t][[1]]]

The resonant solution is the case: ωapp -> ωchar

ResonantSolution =Simplify[yParticularSol/.ωapp -> ωchar]

To analyze the resonant behavior when the viscosity η is small, the solution is expanded about small viscosity to second order:
(Below, the function Map is used with a pure function to simplify each of the terms resulting from a series expansion (n.b. Normal converts the series form into a normal form).   Each term in the expansion is operated on in the following manner:
1) Convert any Exp[] into trignometric functions with ExpToTrig
2) Remove  fractional powers of containing powers of real numbers with PowerExpand
3) Simplify
by using a pure function as the first argument to Map,

ResonantSolutionSmallViscosity =Map[Simplify[PowerExpand [ ExpToTrig[#]]]&,Normal[Series[ResonantSolution,{η,0,2}]]]

ResonantSolutionSmallViscosityDetuned =Map[Simplify[PowerExpand [ ExpToTrig[#]]]&,Normal[Series[yParticularSol,{ωapp,ωchar,1},{η,0,2}]]]

The solution above shows that in the limit of  small viscosity, η,  amplitudes tend to grow linearly with time (i.e.((  Fapp t)/(2 M ωchar)));  the effect of small but finite viscosity η produces "de-tuning" by introducing a time lag.

Of course, this leading behavior could have been obtained directly, viz

ResonatSolZeroViscosity = Simplify[y[t]/.DSolve[{M y''[t] +  Kspring y[t] == Fapp Cos[ωchar t], y[0]==0, y'[0]==0},y[t],t][[1]]]

but, from the direct solution, one cannot obtain the behavior as the viscosity becomes small

Considering the appearance of the resonant solution, there may be another special case when η is near 2 M ωchar.

Map[Simplify[PowerExpand [ #]]&,Normal[Series[ResonantSolution,{η,2 M ωchar,1}]]]

However, it appears that the only special feature is that, after very long times, the solution is a quarter-period out of phase.

Visualizing Solutions

Create a Mathematica  function that returns the solution for specified mass, viscous term, characteristic and applied frequencies

y[M_ , η_ , ωchar_, ωapp_]:= Chop[y[t]/.DSolve[{M y''[t] + η y'[t] + M  ωchar^2 y[t] ==  Cos[ωapp t],y[0]==0,y'[0]==0},y[t],t]//Flatten]

Experiment by plotting for many different values:

Undamped Resonance:

Plot[Evaluate[y[1,0,1/2,1/2]],{t,0,200},PlotPoints->200]

Undamped Near Resonance:

Plot[Evaluate[y[1,0,1/2 +  0.05,1/2]],{t,0,200},PlotPoints->200]

Damped Resonance:

Plot[Evaluate[y[1,1/10,1/2,1/2]],{t,0,200}]

Overdamped Resonance:

Plot[Evaluate[y[1,10,1/2,1/2]],{t,0,200}]

Damped Near Resonance:

Plot[Evaluate[y[1,.05,1/2 +  0.05,1/2]],{t,0,200},PlotPoints->200]

Heavily damped Near Resonance:

Plot[Evaluate[y[1,2.5,1/2 +  0.05,1/2]],{t,0,200},PlotPoints->200]