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  =  
            Random[Real, {-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] +  
                        Random[Real, {-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  =  
      DSolve[M y''[t]  +  η y'[t]  +  
              Kspring y[t]  ==  Fapp Cos[ωapp t], y[t], t] // 
        Flatten

The general solution of the heterogeneous ODE is the sum of the homogeneous solution  (i.e., the one for which the right-hand-side is zero) and the particular solution (i.e., the one for the particular right-hand-side of the ODE)

yParticular[
        0]  =  (y[t] /. yGeneralSol) /. {C[1] -> 0, 
          C[2] -> 0}

yParticular[1]  =  
      Collect[FullSimplify[yParticular[0]], {Sin[t ωapp],  
          Cos[t ωapp]}]

The particular solutions only picks up modes from the forcing term

The homogeneous solution only has he natural frequencies

yHomogenous  = (y[t] /. yGeneralSol)  -  
        yParticular[0]

The General Solution is the combination of two different frequencies--this should give rise to beats if the driving frequency differs from the natural frequency

The following shows that the solution is unphysical when ωapp -> ωchar AND η->0

singbehav  = 
      Series[(yParticular[
                1] /. {η  ->   δ η0 ,  ωapp  ->  ωchar  +  ϵ ω0}), {δ, 0, 1}, {ϵ, 0, 1}] // Normal

We can determine the behavir near the singularity by looking at the effect of each term above:
	1. If the forcing frequency approaches the resonance frequency (ϵ->0), the solution is unbounded because there is no δ in the first part of the second term's numerator
	2. If the viscosity approaches zero (δ->0), the solution will have a linearly growing applitude because of the third term
	3. Other cases will depend on how the ratio δ/ϵ scales.

The following demonstrates the zero viscosity case, needs to be analyzed carefully (as above) otherwise one might miss terms:

yparticularUndamped  =  yParticular[1] /. η -> 0

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] == 1, 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]