MAT 275 Laboratory 6
Forced Equations and Resonance
Exercise 1
Part A
As seen in the graph, the period of forced oscillation for the system is about 5 seconds. The amplitude
was calculated from the equation in the pdf file.
alpha= atan(c*omega/(omega0^2-omega^2)); %alpha derived from the formula
alpha =
0.6015
Part B
function LAB06ex1
clc
omega0 = 2; c = 1; omega = 1.4;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 50;
options = odeset('AbsTol',1e-10,'relTol',1e-10);
[t,Y] = ode45(@f,[t0,tf],Y0,options,param);
y = Y(:,1); v = Y(:,2);
figure(1)
plot(t,y,'b-'); ylabel('y'); grid on;
t1 = 25; i = find(t>t1);
C = (max(Y(i,1))-min(Y(i,1)))/2;
disp(['computed amplitude of forced oscillation = ' num2str(C)]);
Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
disp(['theoretical amplitude = ' num2str(Ctheory)]);
alpha= atan(c*omega/(omega0^2-omega^2)); %alpha derived from the formula
yp=Ctheory*cos(omega*t-alpha); %Ctheory for the equation
yc=y-yp; %yp-particular solution, yc- characteristic equation
figure(2) %creating a new graph
plot(t,yc); grid on %plotting yc versus time
%----------------------------------------------------------------
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 ];
The graph does show exponentially decreasing oscillations because the graph goes towards zero as time
goes towards infinity. The transient function affects the beginning of the graph, but eventually it wears
off and the graph decreases towards zero.
Exercise 2
Part A
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 = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2); %from the equation
end
figure(2)
plot(OMEGA,C, 'ro-', OMEGA, Ctheory); grid on; %plotting omega vs. C and
ctheory
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 ];
The omega value that yields the highest amplitude of the forced oscillation, C, is approximately ω =1.85.
The C value at the practical resonance frequency is approxamtely C= .515.
Part C
C=1/sqrt(OMEGA0^2-OMEGA^2)^2+C^2*OMEGA^2 %the original equation
OMEGA0=2 C=1 %defining the variables
dC/dt=(7*OMEGA-2OMEGA^3)/(OMEGA^4-7*OMEGA^2+16)^(3/2) %the derivative of C
0= (7*OMEGA-2OMEGA^3)/(OMEGA^4-7*OMEGA^2+16)^(3/2) %set the equation = 0
OMEGA= sqrt(7/2)=1.87 %the answer
The OMEGA values that yield the maximum amplitude is 1.87 for the analytically derived formula and
1.85 for the graph. These values are almost the exact same.
Part D
function LAB06ex1
clc
omega0 = 2; c = 0; omega = 2;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 100;
options = odeset('AbsTol',1e-10,'relTol',1e-10);
[t,Y] = ode45(@f,[t0,tf],Y0,options,param);
y = Y(:,1); v = Y(:,2);
figure(1)
plot(t,y,'b-'); ylabel('y'); grid on; hold on
t1 = 25; i = find(t>t1);
C = (max(Y(i,1))-min(Y(i,1)))/2;
disp(['computed amplitude of forced oscillation = ' num2str(C)]);
Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
disp(['theoretical amplitude = ' num2str(Ctheory)]);
alpha= atan(c*omega/(omega0^2-omega^2)) %alpha derived from the formula
yp=Ctheory*cos(omega*t-alpha); %Ctheory for the equation
yc=y-yp; %yp-particular solution, yc- characteristic equation
figure(2) %creating a new graph
plot(t,yc); grid on %plotting yc versus time
%----------------------------------------------------------------
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 ];
This graph has amplitude=1.7 for ω =2. The first graph had amplitude= 0.4 for ω =1.4. The amplitude
increased significantly for the change in ω. If you run LAB06ex1.m with any other value of ω, a direct
correlation to the changing amplitude could be expected for a change in ω.
Part E
C=1 C=2
C=0
OMEGA0 = 2 OMEGA0=1
Changing the OMEGA0 values and the C values in the LAB06ex2.m file changes the results, changing the
OMEGA value only changes the number of points the curve goes through. This makes sense because the
OMEGA0 values and C values are the constants in the file, with a changing OMEGA value, so the OMEGA
is dependent variable that is affected by a change in OMEGA0, C, and time.
Exercise 3
Part A
function LAB06ex2
omega0 = 2; c = 0;
OMEGA = 1:0.02:10;
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 = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2); %from the equation
end
figure(2)
plot(OMEGA,C, 'ro-', OMEGA, Ctheory); grid on; %plotting omega vs. C and
ctheory
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 ];
At ω=2, the spring hits its resonant frequency, causing it to oscillate with a greater amplitude each
oscillation. The maximum amplitude is C=12 , and the OMEGA value and the omega0 value ar the same,
2.
Part B
function LAB06ex1
clc
omega0 = 1; c = 1; omega = 2;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 50;
options = odeset('AbsTol',1e-10,'relTol',1e-10);
[t,Y] = ode45(@f,[t0,tf],Y0,options,param);
y = Y(:,1); v = Y(:,2);
figure(1)
plot(t,y,'b-'); ylabel('y'); grid on;
t1 = 25; i = find(t>t1);
C = (max(Y(i,1))-min(Y(i,1)))/2;
disp(['computed amplitude of forced oscillation = ' num2str(C)]);
Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
disp(['theoretical amplitude = ' num2str(Ctheory)]);
alpha= atan(c*omega/(omega0^2-omega^2)) %alpha derived from the formula
yp=Ctheory*cos(omega*t-alpha); %Ctheory for the equation
yc=y-yp; %yp-particular solution, yc- characteristic equation
figure(2) %creating a new graph
plot(t,yc); grid on %plotting yc versus time
%----------------------------------------------------------------
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 ];
Y vs. t
Plot- characteristic equation
The system stays at maximum amplitude the entire graph, because the spring is as its resonance
frequency, due to the fact that ω=2. The spring cannot oscillate with more amplitude.
Exercise 3
Part A
function LAB06ex1
clc
omega0 = 2; c = 0; omega = 1.8;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 100;
options = odeset('AbsTol',1e-10,'relTol',1e-10);
[t,Y] = ode45(@f,[t0,tf],Y0,options,param);
y = Y(:,1); v = Y(:,2);
figure(1)
plot(t,y,'b-'); ylabel('y'); grid on; hold on
%t1 = 25; i = find(t>t1);
%C = (max(Y(i,1))-min(Y(i,1)))/2;
%disp(['computed amplitude of forced oscillation = ' num2str(C)]);
%Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
%disp(['theoretical amplitude = ' num2str(Ctheory)]);
A=2/abs(omega0^2-omega^2)*sin(.5*(omega0-omega)*t)
plot(t,A,'r',t,-A,'g')
alpha= atan(c*omega/(omega0^2-omega^2)) %alpha derived from the formula
yp=Ctheory*cos(omega*t-alpha); %Ctheory for the equation
yc=y-yp; %yp-particular solution, yc- characteristic equation
figure(2) %creating a new graph
plot(t,yc); grid on %plotting yc versus time
%----------------------------------------------------------------
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 ];
Part B
The period of oscillation for the fast curve is 3.3, as seen in the graph below.
Part C
As seen in the graph it take about 32 seconds to complete a beat. The analytic formula expresses the
time it take to complete half of a beat, so it makes sense that the analytic yields 15.7 seconds.
2*C*sin(.5(OMEGA0OMEGA)t) %the envelope function
OMEGA0=2 OMEGA=1.8 %the values given
C=1/abs(OMEGA0-OMEGA) %the definition of C
y(t)=50/19*sin(.1*t) %plugging in the values given
dy/dt=10*cos(.1*t) %the derivative of the envelope function
0=5/19*cos(.1t) %set it equal to 0
t=15.7 %solve for t
Part D
ω =1.9
ω=1.6
When ω=1.9, the period of forced oscillation is about 63. When ω=1.6, the period of forced oscillation is
about 15. As ω approaches
ω0
, the beat period gets greater, because the system approached the
resonance frequency. The further away from
ω0
that ω gets, the shorter the period of oscillation
gets.
Part E
ω =0.5
The beats phenomenon no longer occurs because the ω is too far away from
ω0
. The ω is not
approaching the resonance frequency at the current value, so the oscillations are happening too fast to
be able to create the beats phenomenon.
