Gain and Phase Shift

Part 1. Background: A Damped Forced Oscillator

Explain the independence of M and ϕ from y[0] and (dy/dt) at 0.

Answer:

Explain the formula for ϕ.

Answer:

Part 2. Finding Gain and Phase Shift

The following lines define an initial value problem and plot the solution (blue) along with the driving function (red).  You may need to enlarge the graph.  To find the coordinates of specific points, click on the graph to select it.  Then, hold down Ctrl and move the cursor onto the graph.  The coordinates will be displayed on the status bar in the lower-left corner of the notebook window.  For more information go to Input->Get Graphics Coordinates in the toolbar above.

omega  =  0.5;
de1  =  y''[t] + 
            Sqrt[1156]*y'[t]/100 + (23/4)^3*y[t]/100 == 
          Cos[omega*t];
solution  =  
        NDSolve[{de1,  y[0] == 0, y'[0] == 0}, 
          y, {t, 0, 40}];
solutionPlot  =  
        Plot[Evaluate[y[t] /. solution], {t, 0, 40}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
input  =  
        Plot[Cos[omega*t], {t, 0, 40}, PlotStyle -> RGBColor[1, 0, 0],  
          DisplayFunction -> Identity];
Show[solutionPlot, input, DisplayFunction -> $DisplayFunction, 
        AxesLabel -> {"t", "y[t]"}];

Record your observations of M[0.5] and ϕ[0.5] in the cell below.

Change the initial conditions and recompute.  Note that we first write a modified ArcCot that uses [0, π] as its range instead of the principle branch.

mArcCot[x_] := If[x < 0, ArcCot[x] + Pi, ArcCot[x]];

c  =  17/100;
k  =  N[Sqrt[12167/6400]];
M[omega_] := 
        1/Sqrt[(k^2 - omega^2)^2 + 4*c^2*omega^2];
phi[omega_] :=  
        mArcCot[(k^2 - omega^2)/(2*c*omega)];

Now change the value of ω successively to 0.25, 0.75, etc., and find M and ϕ for each ω.  A vector of values of ω is already filled in.  Determine the commands necessary to calculate Mdata and phiData.

omegaData  =  {0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0};
Mdata  =  
        Table[M[omegaData[[i]]], {i, Length[omegaData]}];
phiData  =  
        Table[phi[omegaData[[i]]], {i, Length[omegaData]}];

The following commands will plot your data points for M and ϕ.

MplotData  =  Thread[{omegaData, Mdata}]
ListPlot[MplotData, AxesLabel -> {"omega", "M"}, 
        PlotStyle -> {RGBColor[0, 0, 1], 
            PointSize[0.015]}];

phiPlotData  =  Thread[{omegaData, phiData}]
ListPlot[phiPlotData, AxesLabel -> {"omega", "phi"}, 
        PlotStyle -> {RGBColor[0, 0, 1], 
            PointSize[0.015]}];

Part 3. Parameter Identification

Enter the following lines to define the functions M and ϕ.

Clear[k, c];
M[omega_] :=  
        1/Sqrt[(k^2 - omega^2)^2 + 4*c^2*omega^2];
phi[omega_] :=  
        ArcCot[(k^2 - omega^2)/(2*c*omega)];

Decide how to find c and k, and carry out your plan.  You may find the Solve command useful.  The syntax is Solve[{equations,...}, {vars,...}].  We have left ϕ in terms of ArcCot instead of mArcCot so that Solve will recognize it.

Enter your values for c and k, and graph the functions.

c  =  ?;
k  =  ?;
M[omega_] :=  
        1/Sqrt[(k^2 - omega^2)^2 + 4*c^2*omega^2];
phi[omega_] :=  
        mArcCot[(k^2 - omega^2)/(2*c*omega)];

Plot[M[omega], {omega, 0, 3},  
        AxesLabel -> {"omega", "M"}];

Plot[phi[omega], {omega, 0, 3}, 
        AxesLabel -> {"omega", "phi"}];

Part 4. Resonant Frequency

Here again is the definition of M[ω] in terms of parameters k and c.  Decide how locate the maximum value of M[ω].

Clear[k, c, omega];
M[omega_] :=  
        1/Sqrt[(k^2 - omega^2)^2 + 4*c^2*omega^2];

Answer to question 3:

Part 5. Summary

Answers to the questions in Part 5: