Matlab assignment help

MATLAB sessions: Laboratory 6

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 + c

dy

dt + ω20y = cosωt. (L6.1)

This equation models the movement of a mass-spring system similar to the one described in Laboratory 5. The forcing term on the right-hand side of (L6.1) models a vibration, with amplitude 1 and frequency ω (in radians per second = 12π 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 =

√ ω20 − c2/4. When c < 2ω0 the general solution of (L6.1) is

y(t) = e− 1 2 ct(c1 cos(ω1t) + c2 sin(ω1t)) + C cos (ωt− α) (L6.2)

with

C = 1√

(ω20 − ω2) 2 + c2ω2

, (L6.3)

α =

 arctan (

cω ω20−ω2

) if ω0 > ω

π + arctan (

cω ω20−ω2

) if ω0 < ω

(L6.4)

and c1 and c2 determined by the initial conditions. The first term in (L6.2) represents the complementary 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 t due to the decreasing exponential factor. The solution then settles into a (forced) oscillation with amplitude C given by (L6.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 = 2 and c = 1 are fixed. Initial conditions are set to 0. For each value of ω, the amplitude C can 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 ω = 1.4).

1 function LAB06ex1

2 omega0 = 2; c = 1; omega = 1.4;

3 param = [omega0,c,omega];

4 t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 50;

5 options = odeset(’AbsTol’,1e-10,’RelTol’,1e-10);

6 [t,Y] = ode45(@f,[t0,tf],Y0,options,param);

7 y = Y(:,1); v = Y(:,2);

8 figure(1)

9 plot(t,y,’b-’); ylabel(’y’); grid on;

c⃝2011 Stefania Tracogna, SoMSS, ASU 1

MATLAB sessions: Laboratory 6

10 t1 = 25; i = find(t>t1);

11 C = (max(Y(i,1))-min(Y(i,1)))/2;

12 disp([’computed amplitude of forced oscillation = ’ num2str(C)]);

13 Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);

14 disp([’theoretical 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 ];

When executing this program we get

computed amplitude of forced oscillation = 0.40417

theoretical amplitude = 0.40417

Lines 10-14 deserve some explanation. Line 10 defines a time t1 after which we think the contribution of the first term in (L6.2) has become negligible compared to the second term. This depends of course

on the parameter values, in particular c. With c = 1 we obtain e− 1 2 ct ≃ 3.7 × 10−6 for t = 25, so this

is certainly small enough compared to the amplitude seen on Figure L6a. The index i of time values larger than t1 is then determined. The quantity Y(i,1) refers to the values of y associated to times larger than t1 only. The computed amplitude is simply half the difference between the max and the min values. This value is compared to the theoretical value (L6.3).

0 10 20 30 40 50 −0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

y

Figure L6a: Forced oscillation.

1. (a) What is the period of the forced oscillation? What is the numerical value (modulo 2π) of the angle α defined by (L6.4)?

(b) In this question you are asked to modify the file LAB06ex1.m in order to plot the complemen- tary solution of (L6.1), that is, the first term in (L6.2). First define in the file the angle α (alpha) using (L6.4), then evaluate the complementary solution yc by subtracting the quan- tity C cos(ω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 C as a function of ω. We use again ω0 = 2, c = 1 and y(0) = y ′(0) = 0. The

previous problem determined C for a specific value of ω. Here we consider a range of values for ω

c⃝2011 Stefania Tracogna, SoMSS, ASU 2

MATLAB sessions: Laboratory 6

and determine numerically the corresponding amplitude C. We then plot the result as a function of ω, together with the theoretical amplitude from (L6.3). You may need the following MATLAB program.

function LAB06ex2

omega0 = 2; c = 1;

OMEGA = 1:0.02:3;

C = zeros(size(OMEGA));

Ctheory = zeros(size(OMEGA));

t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 50; t1 = 25;

for k = 1:length(OMEGA)

omega = OMEGA(k);

param = [omega0,c,omega];

[t,Y] = ode45(@f,[t0,tf],Y0,[],param);

i = find(t>t1);

C(k) = (max(Y(i,1))-min(Y(i,1)))/2;

Ctheory(k) = ??; % FILL-IN

end

figure(2)

plot(??); grid on; % FILL-IN

xlabel(’\omega’); ylabel(’C’);

%———————————————————

function dYdt = f(t,Y,param)

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 ];

1 1.5 2 2.5 3

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

ω

C

Figure L6b: Amplitude as a function of ω

(a) Fill in the missing parts in the M-file LAB06ex2.m and execute it. You should get a figure like Figure L6b. Include the modified M-file in your lab report.

(b) Examine the graph obtained by running LAB06ex2.m and determine for what (approximate) value of ω the amplitude of the forced oscillation, C, is maximal. This value of ω is called the practical resonance frequency. Give the corresponding maximum value of C.

c⃝2011 Stefania Tracogna, SoMSS, ASU 3

MATLAB sessions: Laboratory 6

(c) Determine analytically the value of ω for which the amplitude of the forced oscillation, C, is maximal by differentiating the expression for C in (L6.3) as a function of ω. Compare the value you find with the value obtained in part (b).

(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 nu- merically (by modifying the initial conditions in LAB06ex2.m) and theoretically (by analyzing the expression for C in (L6.3)). Note that the initial conditions for the DE are y0 and v0.

Resonance

We now investigate what happens to the solution (L6.2), and more specifically to the maximal amplitude C of 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.

1 1.5 2 2.5 3 0

2

4

6

8

10

12

14

ω

C

computed numerically theoretical

(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.

Beats

When c = 0 and ω ̸= ω0 , the solution (L6.2) to (L6.1) reduces to

y(t) = c1 cos(ω0t) + c2 sin(ω0t) + C cos(ωt− α)

with C = 1|ω20−ω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) = 2C sin

( 1

2 (ω0 − ω)t

) sin

( 1

2 (ω0 + ω)t

) .

c⃝2011 Stefania Tracogna, SoMSS, ASU 4

MATLAB sessions: Laboratory 6

When ω is close to ω0 we have that ω + ω0 is large in comparison to |ω0 − ω|. Then sin ( 1 2 (ω0 + ω)t

) is a very rapidly varying function, whereas sin

( 1 2 (ω0 − ω)t

) is a slowly varying function. If we define

A(t) = 2C sin ( 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 ω = 1.8 in LAB06ex1. Also extend the interval of simulation to 100.

0 20 40 60 80 100 −3

−2

−1

0

1

2

3

y

(a) In LAB06ex1 define the “envelope” function A = 2C sin ( 1 2 (ω0 − ω)t

) with C = 1|ω20−ω2|

.

Plot A in red and −A in green, together with the solution. You should obtain Figure L6c. Include the modified M-file.

0 20 40 60 80 100 −3

−2

−1

0

1

2

3

y

Figure L6c: Solution and envelope functions

c⃝2011 Stefania Tracogna, SoMSS, ASU 5

MATLAB sessions: Laboratory 6

(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. Include a graph to support your answer.

(c) What is the length of the beats? Determine the length analytically using the envelope func- tions, and numerically from the graph.

(d) Change the value of ω in LAB06ex1 to 1.9 (a value closer to ω0) and then ω = 1.6 ( 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 ω = 0.5, is the beats phenomenon still present? Why or why not?