MAT 275 Laboratory 6
Forced Equations and Resonance
In this laboratory we take a deeper look at second-order nonhomogeneous equations. We will
concentrate on equations with a periodic harmonic forcing term. This will lead to a study of
the phenomenon known as resonance. The equation we consider has the form
d2y
dt2+cdy
dt +ω2
0y= cos ωt. (1)
This equation models the motion of a mass-spring system similar to the one described in
Laboratory 5. The forcing term on the right-hand side of (1) models a vibration, with amplitude
1 and frequency ω(in radians per second = 1
2πrotation per second = 60
2πrotations per minute,
or RPM) of the plate holding the mass-spring system. All physical constants are assumed to
be positive.
Let ω1=pω2
0−c2/4. When c < 2ω0the general solution of (1) is
y(t) = e−1
2ct(c1cos(ω1t) + c2sin(ω1t)) + Ccos (ωt −α) (2)
with
C=1
qω2
0−ω22+c2ω2
,(3)
α=
arctan cω
ω2
0−ω2if ω0> ω
π+ arctan cω
ω2
0−ω2if ω0< ω (4)
and c1and c2determined by the initial conditions. The first term in (2) represents the com-
plementary solution, that is, the general solution to the homogeneous equation (independent of
ω), while the second term represents a particular solution of the full ODE.
Note that when c > 0 the first term vanishes for large tdue to the decreasing exponential
factor. The solution then settles into a (forced) oscillation with amplitude Cgiven by (3). The
objectives of this laboratory are then to understand
1. the effect of the forcing term on the behavior of the solution for different values of ω, in
particular on the amplitude of the solution.
2. the phenomena of resonance and beats in the absence of friction.
The Amplitude of Forced Oscillations
We assume here that ω0= 4 and c= 2 are fixed. Initial conditions are set to 0. For each value
of ω, the amplitude Ccan be obtained numerically by taking half the difference between the
highs and the lows of the solution computed with a MATLAB ODE solver after a sufficiently
large time, as follows: (note that in the M-file below we set ω= 3.5).
1
LAB06ex1.m
1clear all ; % this deletes all varia bl es
2omega0 = 4; c = 2; omega = 3.5;
3param = [ omega0 ,c , omega ];
4t0 = 0; y0 = 0; v0 = 0; Y0 = [y0 ;v0 ]; tf = 40;
5optio ns = o de set ( 'AbsTol',1e -10 , 'relTol',1e -10);
6[t , Y] = ode45 (@ f ,[ t0 , tf ] , Y0 , op tions , param );
7y = Y (: ,1); v = Y(: ,2);
8figure (1)
9plot (t ,y , 'b - '); ylabel ( 'y'); gr id on ;
10 t1 = 13; i = find ( t > t1 );
11 C = ( max ( Y(i ,1)) - min ( Y(i , 1) ))/2;
12 disp ([ 'compu ted amplitude of forced osci lla tio n = ', num2str (C )]);
13 Ctheory = 1/ sqrt (( omega0 ^2 - omega ^2)^2+( c * omega )^2);
14 disp ([ 'theore tic al amplitude = ', num2str ( Ctheory )]);
15 % -- --- --- --- --- --- --- - -- - -- - -- - -- - -- - -- - -- - --- --- --- --- --- --- --- -
16 function dYdt = f (t , Y , param )
17 y = Y (1); v = Y (2);
18 omega0 = param (1); c = param (2); omega = param (3);
19 dYdt = [ v ; cos ( omega *t) - omega0 ^2* y - c*v ];
20 end
When executing this program we get the plot in Figure 1. and the following output in the
command window:
computed amplitude of forced oscillation = 0.12593
theoretical amplitude = 0.12593
Figure 1: Forced oscillation.
Lines 10-14 deserve some explanation. Line 10 defines a time t1 after which we think the
contribution of the first term in (2) has become negligible compared to the second term. This
depends of course on the parameter values, in particular c. With c= 2 we obtain e−1
2ct ≈
2.3×10−6for t= 13, so this is certainly small enough compared to the amplitude seen on
Figure 1. The index iof time values larger than t1 is then determined. The quantity Y(i,1)
refers to the values of yassociated to times larger than t1 only. The computed amplitude is
2
simply half the difference between the max and the min values. This value is compared to the
theoretical value (3).
1. (a) What is the period of the forced oscillation? What is the numerical value (modulo
2π) of the angle αdefined by (4)?
(b) In this question you are asked to modify the file LAB06ex1.m in order to plot the
complementary solution of (1), that is, the first term in (2). First define in the
file the angle α(alpha) using (4), then evaluate the complementary solution yc by
subtracting the quantity Ccos(ωt −α) from the numerical solution y. Plot the
resulting quantity. Does it look like an exponentially decreasing oscillation? Why or
why not? Include the modified M-file and the corresponding plot.
2. We now consider Cas a function of ω. We use again ω0= 4, c = 2 and y(0) = y0(0) = 0.
The previous problem determined Cfor a specific value of ω. Here we consider a range
of values for ωand determine numerically the corresponding amplitude C. We then plot
the result as a function of ω, together with the theoretical amplitude from (3). The plot
is shown in Figure 2. You may need the following MATLAB program.
LAB06ex2
clear all ; % this delet es all v ariables
omega0 = 4; c = 2;
OMEGA =2:0.01 : 5 ;
C = zeros ( size ( OMEGA ));
Cthe ory = zeros ( size ( OMEGA ));
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0 ; v0 ]; tf =40; t1 = 13;
for k = 1: length ( OMEGA )
omega = OMEGA ( k );
par am = [ ome ga 0 , c , ome ga ];
[t ,Y ] = od e4 5 (@f ,[ t0 , tf ] , Y0 ,[] , p ar am );
i = f ind (t > t1 );
C( k) = ( max ( Y(i ,1)) - min ( Y(i , 1)))/2;
Ctheory (k ) = ??; % FILL - IN the f ormula for Ctheory
end
figure (2)
pl ot (??) ; grid on ; % FILL - IN to plo t C a nd Ct heo ry as a fu nc ti on of OM EG A
xlabel('\ omega '); yl ab el ( 'C');
legend('compu ted numerically ','t h e o r e t ical ')
% --- - --- - --- - --- - --- - --- - --- - --- - --- - --- - --- - -- - - -- - --- - --
func tion dYd t = f (t ,Y , p aram )
y = Y (1); v = Y (2);
omega0 = param (1); c = param (2); omega = param (3);
dYdt = [ v ; cos ( omega *t)- omega0 ^2* y -c * v ];
end
(a) Fill in the missing parts in the M-file LAB06ex2.m and execute it. You should get a
figure like Figure 2. Include the modified M-file in your lab report.
(b) Examine the graph obtained by running LAB06ex2.m and determine for what (ap-
proximate) value of ωthe amplitude of the forced oscillation, C, is maximal. This
value of ωis called the practical resonance frequency. Give the corresponding maxi-
mum value of C.
(c) Determine analytically the value of ωfor which the amplitude of the forced oscilla-
tion, C, is maximal by differentiating the expression for Cin (3) as a function of ω.
Compare this value of ωwith the value obtained in part (b).
3
Figure 2: Amplitude as a function of ω
(d) Run LAB06ex1.m with the value of ωfound in part (c) (include the graph). What is
the amplitude of the forced oscillation? How does it compare with the amplitude of
the forced oscillation in problem 1? If you run LAB06ex1.m with any other value of
ω, how do you expect the amplitude of the solution to be?
(e) Are the results affected by changes in the initial conditions? Answer this question
both numerically (by modifying the initial conditions in LAB06ex2.m) and theoreti-
cally (by analyzing the expression for Cin (3)). Note that the initial conditions for
the DE are y0and v0.
Resonance
We now investigate what happens to the solution (2), and more specifically to the maximal
amplitude Cof the forced oscillation, when we let c→0. The value of ωcorresponding to this
maximal amplitude is called pure resonance frequency. When a mechanical system is stimulated
by an external force operating at this frequency the system is said to be resonant.
3. Set c= 0 in LAB06ex2.m.
(a) Explain what happens. What is the maximal amplitude? What is the value of ω
yielding the maximal amplitude in the forced solution? How does this value compare
to ω0?
(b) Run LAB06ex1.m with c= 0 and ωequal to the value found in part (a). Comment
on the behavior of the solution. Include the graph.
4
Beats
When c= 0 and ω6=ω0, the solution (2) to (1) reduces to
y(t) = c1cos(ω0t) + c2sin(ω0t) + Ccos(ωt −α)
with C=1
|ω2
0−ω2|. If the initial conditions are set to zero, the solution reduces to
y(t) = C(cos(ωt)−cos(ω0t))
which can be rewritten as
y(t)=2Csin 1
2(ω0−ω)tsin 1
2(ω0+ω)t.
When ωis close to ω0we have that ω+ω0is large in comparison to |ω0−ω|. Then sin 1
2(ω0+ω)t
is a very rapidly varying function, whereas sin 1
2(ω0−ω)tis a slowly varying function. If we
define A(t) = 2Csin 1
2(ω0−ω)t, then the solution can be written as
y(t) = A(t) sin 1
2(ω0+ω)t
and we may interpret it as a rapidly oscillating function with period T=4π
ω0+ω, but with a
slowly varying amplitude A(t). This is the phenomenon known as beats. Note that A(t) and
−A(t) are the so-called “envelope functions”. The period of A(t) is 4π
|ω0−ω|, thus the length of
the beats is 2π
|ω0−ω|.
4. To see the beats phenomenon, set c= 0 and ω= 3.8 in LAB06ex1. Also extend the
interval of simulation to 80.
(a) In LAB06ex1 define the “envelope” function A= 2Csin 1
2(ω0−ω)twith C=
1
|ω2
0−ω2|.
Plot Ain red and −Ain green, together with the solution. You should obtain Figure
3. Include the modified M-file.
(b) What is the period of the fast oscillation (that is, the period of sin 1
2(ω0+ω)t)?
Confirm your answer by zooming in on the graph of the solution.
(c) What is the length of the beats? Determine the length analytically using the envelope
functions, and numerically from the graph.
5
Figure 3: Solution and envelope functions
(d) Change the value of ωin LAB06ex1 to 3.9 (a value closer to ω0) and then ω= 3.5 (a
value farther away from ω0). Include the two graphs. For each of these two values
of ωfind the period of the fast oscillation and the length of the beats. How do the
periods change compared to parts (b) and (c)?
(e) If you let ω= 2, is the beats phenomenon still present? Why or why not?
6
