Project Help Notebook by Keith Stroyan

Help NoteBook Copyright 1998 by Academic Press, Inc. - All rights reserved.
Modify this Notebook to prepare your project report.

Project 37: Resonance, The Big Wiggle

by A. Student and B. Pupil

with technical help from Da LabMonitor and Help Notebook

READ THE PRINT PROJECT.  This notebook is only to help you understand that and reinforce it with Mathematica computaitons.

Initialization

Needs["Graphics`Colors`"];

My Project Introduction

Write this last.  Begin with the basic question and state your main finding clearly, but with few technical details.  You might even modify the title of this section to make it more interesting.

Natural Oscillation

The autonomous (also called homogeneous) equation with zero external forcing,

m y" + c y' + s y == 0

has solutions of the form

(Subscript[y, T])[t] =  
      (e^(-α t)) 
        ((Subscript[k, 1]) Cos[ν t] + (Subscript[k, 2 ]) Sin[ν t])

where

ν = ((Sqrt[4 m s - (c^2)])/(2 m))  ("nu" Greek "n" as in natural.)

α = (c/(2 m))

provided the term inside the square root is positive.  THIS IS BECAUSE ... (see Theorem 23.5 ff. of text.)

This is the case when the autonomous equation has oscillatory solutions at frequency  ν.  The frequency ν  is called the "natural frequency" of the system, because it oscillates at this frequency by itself, without external forcing.

GIVE SOME SPECIFIC EXAMPLES HERE ...

Damped Oscillatory Solutions

...

Clear[y, t, α, ν, k]; 
    (Subscript[y, T])[t_] =  
      (E^(-α t)) 
        ((Subscript[k, 1]) Cos[ν t] + (Subscript[k, 2 ]) Sin[ν t])
ν = ((Sqrt[4 m s - (c^2)])/(2 m));
α = (c/(2 m));
(Subscript[y, T])[t_]

...

Clear[x, t, m, c, s]; m = 1; c = 1/4; s = 2;
Print["The general solution to ", m, "x''[t] + ", c, 
        "x'[t] + ", s, "x[t] = 0"];
DSolve[{m x''[t]+c x'[t]+s x[t]==0},x[t],t]

Typical graphs of the solutions look as follows.

soln=DSolve[{x[0] == 1,x'[0] == 1,m x''[t]+c x'[t]+s x[t]==0},x[t],t];
x[t_] = x[t] /. soln;
Plot[x[t], {t, 0, 10}, PlotStyle -> Red];

THE NATURAL FREQUENCY OF THIS SYSTEM IS:

ν

Over Damped Solutions

Clear[x, t, a, b, c]; m = 1; c = 4; s = 2;
Print["The general solution to ", m, "x''[t] + ", c, 
        "x'[t] + ", s, "x[t] = 0"];
DSolve[{m x''[t]+c x'[t]+s x[t]==0},x[t],t]

Typical graphs of the solutions look as follows.

soln=DSolve[{x[0] == 1,x'[0] == 1,m x''[t]+c x'[t]+s x[t]==0},x[t],t];
x[t_] = x[t] /. soln;
Plot[x[t], {t, 0, 10}, PlotStyle -> Red];

THIS SYSTEM HAS NO NATURAL FREQUENCY BECAUSE...

Transients just pass through
	(lim_(t -> ∞)) (Subscript[y, T])[t] =   0

The solutions to the equation with zero external forcing are called "transient" because they "pass through the system" or tend to zero,

(lim_(t -> ∞)) (Subscript[y, T])[t] =   0, because ...

Sinusoidal Forcing

If we apply a sinusoidal force

F[t] = (Subscript[f, C]) Cos[ω t] +  (Subscript[f, S]) Sin[ω t]

to a linear equation, say a spring - shock absorber - mass system,

m y''[t] + c y'[t] + s y[t]=F[ t]

m y''[t] + c y'[t] + s y[t]=(Subscript[f, C]) Cos[ω t]+ (Subscript[f, S]) Sin[ω t]

then a solution with any initial condition tends toward the sinusoidal "stead state" solution,

(Subscript[y, SS])[t]  =  h Cos[ω t]  +  k Sin[ω t]

where the constants  h  and  k  are given by

h = ((det|({{(Subscript[f, C]), c ω}, {(Subscript[f, S]), s  -  m (ω^2)}})|)/(det|({{s  -  m (ω^2), c ω}, {-c ω, s  -  m (ω^2)}})|))  and  k = ((det|({{s  -  m (ω^2), (Subscript[f, C])}, {-c ω, (Subscript[f, S])}})|)/(det|({{s  -  m (ω^2), c ω}, {-c ω, s  -  m (ω^2)}})|))

Why does any solution tend to the steady state oscillation?

(See Project 36 and notice that because  c > 0, the solutions to the "homogenous" or "autonomous" problem with zero applied force consists of sums of exponentials with negative powers.  These tend to zero.  Use Superposition to write any solution as a sum of these basic zero force solutions and the steady state solution.)

Every solution to the forced equation can be wirtten in the form

y[t]  =  (Subscript[y, T])[t] +  (Subscript[y, SS])[t]

where  (Subscript[y, SS])[t]  is the steady state solution and  (Subscript[y, T])[t]  satisfies the un-forced equation  m y''[t]+c y'[t]+s y[t] = 0.

...

Example

In the naturally oscillating case, ...

y[t]  =  (Subscript[y, T])[t] +  (Subscript[y, SS])[t]

(Subscript[y, T])[t] =  
      (e^(-α t)) 
        ((Subscript[k, 1]) Cos[ν t] + (Subscript[k, 2 ]) Sin[ν t])

(Subscript[y, SS])[t]  =  h Cos[ω t]  +  k Sin[ω t]

...

The special case  m y''[t] + c y'[t] + s y[t]=Sin[ω t]

The particular case when the applied force is  Sin[ω t], the steady state solution of the differential equation is

(Subscript[y, SS])[t]  =  h Cos[ω t]  +  k Sin[ω t]

where \(h\  = \ 
        \(c\ ω\)\/\(\((s\  - \ m\ ω\^2)\)\^2 + \((c\ ω)\)\^2\)\)  and  \(k\  = \ 
        \(s\  - \ m\ ω\^2\)\/\(\((s\  - \ m\ ω\^2)\)\^2 + \((c\ ω)\)\^2\)\)

because...

Phase and Amplitude of  (Subscript[y, SS])[t]  =  h Cos[ω t]  +  k Sin[ω t]= a Sin[ω t + ϕ]

Given  h  and  k, we may write the solution in the form

h Cos[ω t] + k Sin[ω t] = a Sin[ω t + ϕ]

where  a = (Sqrt[(h^2) + (k^2)])  and  ϕ  is the angle that makes the unit vector {(h/a), (k/a)}= {Cos[ϕ],Sin[ϕ]}, or  ϕ = ArcTan[k/h], with the correct quadrant.  This is because ....

...

Hint 37.1: The Phase-Amplitude Solution to m y''[t] + c y'[t] + s y[t]=Sin[ω t]

The particular case when the applied force is  Sin[ω t], the steady state solution of the differential equation is

a Sin[ω t + ϕ]

where  a = ???

Clear[m, s, c, ω];
h  =  
        ((c ω)/((((s  -  m (ω^2)))^2) + (((c ω))^2)));
k  =  
        ((s  -  m (ω^2))/((((s  -  m (ω^2)))^2) + (((c ω))^2)));
a = (Sqrt[(h^2) + (k^2)])
Simplify[a]

because...

Amplitude of the Steady State Solution 
as a Function of forcing Frequency

Hint 37.2

The maximum amplitude of the steady state solution,

Max[ a[ω] : ω ≥ 0]

occurs at the same frequency that minimizes

Min[(((s  -  m (ω^2)))^2) + (((c ω))^2): ω ≥ 0 ]

because ...

Conclude the main part of your project with a CAREFUL computation of the minimum of  b[ω]  or the maximizing frequency for  a[ω].  Following is some help that tries to prevent you from making gross logical errors.

Minimize (((s  -  m (ω^2)))^2) + (((c ω))^2)

Clear[m, c, s, ω, b];
b[ω_] = (((s  -  m (ω^2)))^2) + (((c ω))^2)

b'[ω]

Solve[b'[ω] == 0, ω]

Oscillation

Clear[y, t, α, ν, k, m, c, s, a];
(Subscript[y, T])[t_] =  
      (E^(-α t)) 
        ((Subscript[k, 1]) Cos[ν t] + (Subscript[k, 2 ]) Sin[ν t])
ν = ((Sqrt[4 m s - (c^2)])/(2 m));
α = (c/(2 m));
(Subscript[y, T])[t_]
a[ω_] := 
        (Sqrt[(((c^2) (ω^2))/((((c^2) (ω^2) + (((s - m (ω^2)))^2)))^2)) + ((((s - m (ω^2)))^2)/((((c^2) (ω^2) + (((s - m (ω^2)))^2)))^2))]);
ωMAX = ((Sqrt[2 m s - (c^2)])/((Sqrt[2]) m));
m = 1; c = 1/4; s = 2;
ν
(Subscript[y, T])[t]
ωMAX
N[ν]
N[ωMAX]

Plot[a[ω], {ω, 0, 10}, PlotStyle -> Red];

Oscillation?

m = 1; c = 1; s = 1/2;
ν
(Subscript[y, T])[t]
ωMAX = ((Sqrt[-(c^2) + 2 m s])/((Sqrt[2]) m))
N[ν]
N[ωMAX]

Plot[a[ω], {ω, 0, 10}, PlotStyle -> Blue];

Maximum?

m = 1; c = 1; s = 1/3;
ν
(Subscript[y, T])[t]
ωMAX = ((Sqrt[-(c^2) + 2 m s])/((Sqrt[2]) m))
N[ν]
N[ωMAX]

Plot[a[ω], {ω, 0, 10}, PlotStyle -> Blue];

Oscillation?

m = 1; c = 1; s = 1/5;
ν
(Subscript[y, T])[t]
ωMAX = ((Sqrt[-(c^2) + 2 m s])/((Sqrt[2]) m))
N[ν]
N[ωMAX]

Plot[a[ω], {ω, 0, 10}, PlotStyle -> Blue];

Animation of Resonance

Once you have CAREFULLY maximized the amplitude of the steady state response to forcing, you can illustrate the resonant frequency by showing an animation of the response as forcing frequency varies.  Following is some help (with a trick in case you're not paying attention).

Clear[y, t, h, k, m, c, s, ωMAX, a];
h  :=  
        ((c ω)/((((s  -  m (ω^2)))^2) + (((c ω))^2)));
k  :=  
        ((s  -  m (ω^2))/((((s  -  m (ω^2)))^2) + (((c ω))^2)));
ωMAX := ((Sqrt[2 m s - (c^2)])/((Sqrt[2]) m));
a[ω_] := 
        (Sqrt[(((c^2) (ω^2))/((((c^2) (ω^2) + (((s - m (ω^2)))^2)))^2)) + ((((s - m (ω^2)))^2)/((((c^2) (ω^2) + (((s - m (ω^2)))^2)))^2))]);
(Subscript[y, SS])[t_]  :=  h Cos[ω t]  +  k Sin[ω t];

m = 1.0;
c = 0.1;
s = 1.0;
For[ω =  .5, ω < 1.5, 
        ω +=  .1, 
Plot[(Subscript[y, SS])[t], {t, 0, 50}, 
          PlotRange -> {-10, 10}, PlotLabel -> ω]];
ωMAX
a[ωMAX]

m = 1.0; c = 0.5; s = 2.0;
For[ω =  .5, ω < 2, 
        ω +=  .1, 
Plot[(Subscript[y, SS])[t], {t, 0, 50}, 
          PlotRange -> {-2, 2}, PlotLabel -> ω]];
ωMAX
a[ωMAX]

m = 1; c = 1; s = 1/3;
For[ω =  .5, ω < 2, 
        ω +=  .1, 
Plot[(Subscript[y, SS])[t], {t, 0, 50}, 
          PlotRange -> {-2, 2}, PlotLabel -> ω]];
ωMAX
a[ωMAX]