Forced Vibrations

forced.mws

Douglas B. Meade (meade@math.sc.edu)

1 Overview

The complete solution to a second-order linear constant coefficient ODE with a trigonometric forcing function can be difficult to determine by inspection. Particularly difficult is the detection of resonance. Several different characteristics of the solution are identified in this section, and this worksheet. The discussion in the text is valid for any problem of this type (i.e., the analysis is done for general parameters). The same analysis could be repeated using Maple. However, the general formulae are only a small part of the story; understanding the physical meaning of the different quantities that will be discussed is, in my opinion, of more long-lasting importance. Thus, the discussion found in this worksheet will generally work with specific sets of parameters. It is hoped that this presentation, combined with the description in the text, combine to create a more complete picture of this important material.

2 Objectives

• introduce concept of frequency gain, frequency response curve, and resonance frequency

• emphasize physical understanding and computations from first principles over memorization of formulae

• apply and combine Maple techniques for solving second-order linear constant coefficient ODEs to extract information from the solution

3 New Maple Commands

• arctan arctangent function (particularly the two-argument form)

• coeff extract the coefficient of a specific term of an expression • collect collect coefficients of specified terms in an expression

• combine(..., trig) combines a product of trigonometric terms into a sum of trigonometric terms

4 A Basic Problem

> restart;
> with( plots ):
> with( DEtools ):

Consider the second-order linear constant coefficient ODE

> ODE := diff( x(t), t,t ) + 2*diff( x(t), t ) + 5*x(t) = 3*cos(4*t);

ODE := (\frac{∂^{2}}{∂t^{2}} x(t)) + 2 (\frac{∂}{∂t} x(t)) + 5 x(t) = 3 cos(4 t)

with homogeneous initial conditions

> IC1 := x(0)=0;
> IC2 := D(x)(0)=0;

IC1 := x(0) = 0

IC2 := D(x)(0) = 0

The equation of motion (i.e., the solution) to this IVP can be obtained directly using dsolve:

> EOM := dsolve( { ODE, IC1, IC2 }, x(t) );

EOM := x(t) = −\frac{1}{740}(−480 sin(2 t) cos(t)^{6} e^{t} − 2880 sin(2 t) cos(t)^{5} sin(t) e^{t} + 720 sin(2 t) cos(t)^{4} e^{t} + 2880 sin(2 t) cos(t)^{3} sin(t) e^{t} − 492 sin(2 t) cos(t)^{2} e^{t} − 984 sin(2 t) cos(t) sin(t) e^{t} + 126 e^{t} sin(2 t) − 90 cos(2 t) cos(6 t) e^{t} + 15 cos(2 t) sin(6 t) e^{t} + 222 cos(2 t)^{2} e^{t} − 111 cos(2 t) sin(2 t) e^{t} + 126 sin(2 t) − 132 cos(2 t))/e^{t} This solution is \underline{NOT} what we would have expected!

There are trigonometric identities that combine a product of trigonometric terms to form a sum of trigonometric terms. The Maple command that applies these identities is combine (with an optional second argument of trig ). See the on-line help for combine and combine,trig for a full description of the uses of these commands.

   > SOLN := combine( rhs(EOM), trig );

SOLN := \frac{1}{370} \frac{−63 sin(2 t) + 66 cos(2 t) + 48 e^{t} sin(4 t) − 66 e^{t} cos(4 t)}{e^{t}}

Further improvements can be achieved using collect (see ?collect):

   > SOLN := collect( SOLN, exp(t) );

SOLN := \frac{24}{185} sin(4 t) − \frac{33}{185} cos(4 t) + \frac{−\frac{63}{370} sin(2 t) + \frac{33}{185} cos(2 t)}{e^{t}}

4.1 Points to Ponder

• What should the solution look like?

• Why?

• How did Maple come up with this solution?

4.2 Analysis of the Solution

4.2.1 Identify the transient and steady-state solutions

4.2.2 Solution

The transient solution consists of all terms that decay to zero as t increases; the steady-state solution contains all other terms. Another way to look at this problem is that the steady-state solution is the particular solution and the transient solution is the homogeneous solution to the IVP.

   > tSOLN := (-63/370*sin(2*t)+33/185*cos(2*t))/exp(t);

tSOLN := \frac{−\frac{63}{370} sin(2 t) + \frac{33}{185} cos(2 t)}{e^{t}}

> ssSOLN := 24/185*sin(4*t)-33/185*cos(4*t);

ssSOLN := \frac{24}{185} sin(4 t) − \frac{33}{185} cos(4 t)

> plot( [ SOLN, tSOLN, ssSOLN ], t=0..5,
> color=[RED,GREEN,BLUE], style=[POINT,LINE,LINE],
> title=`Full, Transient, and Steady-State Solutions` );

4.2.3 What is the amplitude of the forcing term? 4.2.4 Solution A quick look at the right-hand side of the ODE shows that the amplitude is

> A[f] := 3;

A_{f} := 3

4.2.5 What are the amplitude and phase angle of the steady-state

solution? 4.2.6 Solution Recall that the steady-state solution is

> ssSOLN;

\frac{24}{185} sin(4 t) − \frac{33}{185} cos(4 t)

The amplitude of the solution is the square root of the sum of the squares of the coefficient of each of the two terms:

> A[ss] := sqrt( coeff( ssSOLN, sin(4*t) )^2 + coeff( ssSOLN, cos(4*t)
)^2 );
> evalf( A[ss] );

A_{ss} := \frac{3}{185} \sqrt{185}

.2205643867

The phase angle is also determined from these coefficients.

> phi[ss] := arctan( coeff( ssSOLN, cos(4*t) ), coeff( ssSOLN,
sin(4*t) ) );
> evalf( phi[ss] );

ϕ_{ss} := −arctan(\frac{11}{8})
−.9420000404

4.2.7 Note about use of the arctan function Note the use of the two-argument arctan function. The advantage over the one-argument arctan is that angles can be returned in the appropriate quadrant (e.g., arctan(1) can be either \frac{π}{4} or \frac{5 π}{4} but arctan(-2,-2) is definitely \frac{5 π}{4} . (See ?arctan for additional information.)

Thus, an equivalent form for the steady-state solution is

> ssSOLNa := A[ss] * cos( 4*t - phi[ss] );

ssSOLNa := \frac{3}{185} \sqrt{185} cos(4 t + arctan(\frac{11}{8}))

To verify that this actually is the steady-state solution:

> expand( ssSOLNa, trig );

\frac{192}{185} cos(t)^{4} − \frac{192}{185} cos(t)^{2} + \frac{24}{185} − \frac{264}{185} sin(t) cos(t)^{3} + \frac{132}{185} sin(t) cos(t)

> combine( ", trig );

\frac{24}{185} cos(4 t) − \frac{33}{185} sin(4 t)

   > ssSOLN;

\frac{24}{185} sin(4 t) − \frac{33}{185} cos(4 t)

5 Frequency Gain

Definition 1: The FREQUENCY GAIN (or GAIN FACTOR) is the ratio of the amplitudes of the steady-state solution and the forcing term.

5.1 Example

In the current example, the frequency gain is

   > M := A[ss]/A[f];
   > evalf(M);

M := \frac{1}{185} \sqrt{185}

.07352146221

Note that this ratio is, in some sense, the fraction of the forcing term that that is present in the steady-state solution. In some situations it will be desired to maximize the frequency gain of a system.

The book’s discussion of frequency gain points out that the gain factor depends, in part, on the frequency of the forcing term (and is independent of the amplitude of the forcing term). This suggests considering the same problem, without selecting a specific value for the frequency of the forcing function. That is,

   > ODE2 := diff( x(t), t,t ) + 2*diff( x(t), t ) + 5*x(t) = 3*cos(gamma*t);

ODE2 := (\frac{∂^{2}}{∂t^{2}} x(t)) + 2 (\frac{∂}{∂t} x(t)) + 5 x(t) = 3 cos(γ t)

Note that gamma is a defined constant in Maple (it refers to Euler’s constant, see ?ininame). The full set of constants is

   > constants;

false, γ, ∞, true, Catalan, FAIL, π

Constants can be added or removed by the user. In this case, we can remove gamma as follows:

   > constants := op( { constants } minus { gamma } );

constants := ∞, FAIL, false, π, true, Catalan

5.2 Note about protected names (such as gamma)

Even though gamma has been removed from the list of constants, it is still a protected name in Maple. We can use gamma in our calculations, but cannot assign a value to gamma until it is unprotected. (See ?protect for more information.)

> EOM2 := dsolve( { ODE2, IC1, IC2 }, x(t) );

EOM2 := x(t) = \frac{1}{4}(15 sin(2 t) cos((−2 + γ) t) + 3 cos(2 t) sin((−2 + γ) t) γ^{2}

+ 12 cos(2 t) sin((−2 + γ) t) γ − 3 cos(2 t) cos((−2 + γ) t) γ^{3} − 6 cos(2 t) cos((−2 + γ) t) γ^{2} + 9 cos(2 t) cos((−2 + γ) t) γ − 3 cos(2 t) sin((2 + γ) t) γ^{2} + 12 cos(2 t) sin((2 + γ) t) γ + 3 cos(2 t) cos((2 + γ) t) γ^{3} − 6 cos(2 t) cos((2 + γ) t) γ^{2} + 3 sin(2 t) sin((2 + γ) t) γ^{3} − 6 sin(2 t) sin((2 + γ) t) γ^{2} − 9 sin(2 t) sin((2 + γ) t) γ + 3 sin(2 t) cos((2 + γ) t) γ^{2} − 12 sin(2 t) cos((2 + γ) t) γ + 3 sin(2 t) sin((−2 + γ) t) γ^{3} + 6 sin(2 t) sin((−2 + γ) t) γ^{2} − 9 sin(2 t) sin((−2 + γ) t) γ + 3 sin(2 t) cos((−2 + γ) t) γ^{2} + 12 sin(2 t) cos((−2 + γ) t) γ − 9 cos(2 t) cos((2 + γ) t) γ + 30 sin(2 t) sin((2 + γ) t) + 15 sin(2 t) cos((2 + γ) t) − 30 sin(2 t) sin((−2 + γ) t) + 15 cos(2 t) sin((−2 + γ) t) + 30 cos(2 t) cos((−2 + γ) t)

− 15 cos(2 t) sin((2 + γ) t) + 30 cos(2 t) cos((2 + γ) t) − 72 \frac{(γ^{2} − 5) e^{(−t)} cos(2 t) γ^{2}}{%1}

+ 12 \frac{(γ^{2} − 5) e^{(−t)} cos(2 t) γ^{4}}{%1} + 36 \frac{(γ^{2} + 5) e^{(−t)} sin(2 t) γ^{2}}{%1} − 6 \frac{(γ^{2} + 5) e^{(−t)} sin(2 t) γ^{4}}{%1}

− 150 \frac{(γ^{2} + 5) e^{(−t)} sin(2 t)}{%1} + 300 \frac{(γ^{2} − 5) e^{(−t)} cos(2 t)}{%1})/(%1)

%1 := 25 − 6 γ^{2} + γ^{4}

This is even worse than before! Let’s see if combine(..., trig) can transform this into a form we can use.

> SOLN2 := combine( rhs(EOM2), trig );

SOLN2 := (−6 γ^{2} cos(γ t) + 6 e^{(−t)} cos(2 t) γ^{2} − 3 e^{(−t)} sin(2 t) γ^{2} + 12 γ sin(γ t)

− 15 e^{(−t)} sin(2 t) + 30 cos(γ t) − 30 e^{(−t)} cos(2 t)) /(50 − 12 γ^{2} + 2 γ^{4})

To be most effective, it would be nice to have all the terms involving the exponential solution and the different trigonometric terms grouped together. This can be done using the collect command.

> SOLN2 := collect( SOLN2,
> [ exp(-t), sin(2*t), cos(2*t), sin(gamma*t), cos(gamma*t) ] );

SOLN2 :=

(\frac{(−3 γ^{2} − 15) sin(2 t)}{%1} + \frac{(6 γ^{2} − 30) cos(2 t)}{%1}) e^{(−t)} + 12 \frac{γ sin(γ t)}{%1} + \frac{(−6 γ^{2} + 30) cos(γ t)}{%1}

%1 := 50 − 12 γ^{2} + 2 γ^{4}

As before, the transient and steady-state solutions are easily identified (expressions obtained verbatim, except for adding spaces for readability, by cutand-paste from previous results):

> tSOLN2 := ( (-3*gamma^2-15)/(-12*gamma^2+2*gamma^4+50) * sin(2*t)
> + (6*gamma^2-30)/(-12*gamma^2+2*gamma^4+50) * cos(2*t) ) * exp(-t);

tSOLN2 := (\frac{(−3 γ^{2} − 15) sin(2 t)}{50 − 12 γ^{2} + 2 γ^{4}} + \frac{(6 γ^{2} − 30) cos(2 t)}{50 − 12 γ^{2} + 2 γ^{4}}) e^{(−t)}

> ssSOLN2 := 12*gamma/(-12*gamma^2+2*gamma^4+50) * sin(gamma*t)
> + (-6*gamma^2+30)/(-12*gamma^2+2*gamma^4+50) * cos(gamma*t);

ssSOLN2 := 12 \frac{γ sin(γ t)}{50 − 12 γ^{2} + 2 γ^{4}} + \frac{(−6 γ^{2} + 30) cos(γ t)}{50 − 12 γ^{2} + 2 γ^{4}}

The amplitude of the steady-state solution is

> A[ss2] := sqrt( simplify( coeff( ssSOLN2, sin(gamma*t) )^2
> + coeff( ssSOLN2, cos( gamma*t) )^2 ) );

A_{ss2} := 3 \sqrt{\frac{1}{25 − 6 γ^{2} + γ^{4}}}

> phi[ss2] := arctan( coeff( ssSOLN2, cos(gamma*t) ), coeff( ssSOLN2,
sin(gamma*t) ) );

ϕ_{ss2} := arctan(\frac{−6 γ^{2} + 30}{50 − 12 γ^{2} + 2 γ^{4}} , 12 \frac{γ}{50 − 12 γ^{2} + 2 γ^{4}} )

A quick check of these quantities when γ = 4, produces the same quantities as before

> eval( subs( gamma=4, [ A[ss2], phi[ss2] ] ) );

[\frac{3}{185} \sqrt{185}, −arctan(\frac{11}{8})]

Moreover, the frequency gain for each value of gamma is easily computed to be

> M := A[ss2] / A[f] ;

M := \sqrt{\frac{1}{25 − 6 γ^{2} + γ^{4}}}

6 Frequency Response Curvve

Definition 2: The FREQUENCY RESPONSE CURVE or RESONANCE CURVE for a system is a graph of frequency gain versus the frequency of the forcing term.

6.1 Example

In this example, the frequency response curve is

   > plot( M, gamma=0..10, title=`Frequency Response Curve (M vs.
   gamma)` );

Note that the maximum frequency gain is ˜ 0.248 and occurs for γ ˜ 1.7. The exact location of this maximizer can be found using Calculus. To begin, compute and plot the first derivative of the frequency gain.

   > dM := diff( M, gamma );

dM := −\frac{1}{2} \frac{−12 γ + 4 γ^{3}}{\sqrt{\frac{1}{25 − 6 γ^{2} + γ^{4}}} (25 − 6 γ^{2} + γ^{4})^{2}}

   > dM := simplify( dM, power, symbolic );

dM := −2 \frac{γ (−3 + γ^{2})}{(25 − 6 γ^{2} + γ^{4})^{3/2}}

> plot( dM, gamma=0..10, title=`M'(gamma) vs.  gamma` );

The maximum frequency gain occurs for a (positive) frequency where the derivative has a zero. These points are found to be

> solve( dM=0, { gamma } );

\{γ = 0\}, \{γ = \sqrt{3}\}, \{γ = −\sqrt{3}\}

Since we know gamma is positive, the maximum frequency gain occurs when

> gamma_r := "[2];

gamma r := \{γ = \sqrt{3}\}

Thus, the maximum frequency gain is

> simplify( subs( gamma_r, M ) );

\frac{1}{4}

or 25% of the forcing function’s amplitude.

6.1.1 Cautionary Advice

The ordering of Maple’s response to the solve command may change from one run to the next. Thus, it might be necessary to change the index in the previous command. Just be sure that you get the only positive solution.

6.2 More Points to Ponder

• What is the phase angle when the frequency gain is maximized?

• How does the phase angle vary with the frequency, gamma?

• Describe the small and large frequency behavior of the phase angle

7 Resonance Frequency

Definition 3: The RESONANCE FREQUENCY of a system is the frequency, \frac{γ_{r}}{2 π}, at which the maximum frequency response is obtained.

7.1 Example

For the example considered here, resonance occurs when

> RESONANCE := subs( gamma_r, gamma/(2*Pi) );

RESONANCE := \frac{1}{2} \frac{\sqrt{3}}{π}

8 Concluding Remarks

The discussion in the text concludes with a brief investigation of the dependence of the resonance frequency on the damping coefficient. This is omitted from the worksheet, but you are expected to read and understand this relationship. While I encourage you to use Maple to assist with this analysis, you can do the work by hand computations. Either way, you are responsible for this material on the next exam.

9 Copyright

forced.mws – a Maple worksheet illustrate concepts relating to forced vibrations (for use with section 5.3, p. 253, of Nagle and Saff) Prepared by: Douglas B. Meade (7 March 1996) as part of the Maple Supplement for Nagle/Saff (Addison-Wesley) Updated to Maple V, Release 4, by Douglas B. Meade (27 December, 1996) For more information, please contact the author at one of following addresses: Department of Mathematics, Univ. of South Carolina, Columbia, SC 29208 (803) 777-6183 meade@math.sc.edu (copyright, Douglas B. Meade, 1996)