Exercise 1

% Exercise 1(a)
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-','linewidth',1.45); 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)]);
%---------------------------------------------------------------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];

(a)
computed amplitude of forced oscillation = 0.40417
theoretical amplitude = 0.40417
T=

2 2
=
=4.49
1.4

seconds

 =arctan

c
1.4
=arctan 2
2
2
s
2 1.4
2
0

= 0.6015 radians

(b)
% Exercise 1(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);
alpha=atan(c*omega/(omega0^2-omega^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)]);
figure (2)
yc = y-C*cos(omega*t-alpha);
plot (t,yc); grid on
title ('complementary solution')
%---------------------------------------------------------------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 oscillation appears to be decreasing exponentially because the amplitude is decreasing at a

steady pace until the graph ultimately becomes a straight line at x=0.
Exercise 2
% Exercise 2(a)
function LAB06ex2
clc
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) = 1/sqrt((omega0^2-omega^2)^2+c^2*omega^2);
end
figure(2)
plot(OMEGA,C,'b',OMEGA,Ctheory,'ro-'); grid on;
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 ];

(a)

(b)
Zoom in :

Max (practical resonance frequency) = 1.88

Max Amplitude: 0.5163
(c)

C=

( ) +c
2 2

2
0

1/ 2

=( ( 20 2 ) ( c )2)

Set C( ) = 0

( )
1
2

( (20 2 ) +( c )2)

3
2

( 2 (202) (2 ) 2 c2 )=0

(2 ( 202 )+c 2 )=0

=0
c2
2 ( 202 ) + c2=0 => 2 =20
2

= 20

c2
2

0 =2c =1

1
= 4 =1.87
2
In Part B was 1.88 and C was 0.5163, but in Part C we get 1.87 and C = 0.516398 which is
very close to the same.
(d)
% Exercise 2(d)
function LAB06ex1
clc
omega0 = 2; c = 1; omega = 1.87;
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)]);

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

In work space Matlab we get the answer :

computed amplitude of forced oscillation = 0.5164
theoretical amplitude = 0.5164
The amplitude is 0.40417 in after the in problem one, but when the goes to 1.87 the amplitude
of forced oscillation increased to 0.5164. If you decrease the amplitude the oscillation will
decrease and when you increase , the amplitude of oscillation will increase.
(e)
% Exercise 2 (e)
function LAB06ex2
clc
omega0 = 2; c = 1;
OMEGA = 1:0.02:3;
C = zeros(size(OMEGA));
Ctheory = zeros(size(OMEGA));
t0 = 0; y0 = 10; v0 = 12; 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) = 1/sqrt((omega0^2-omega^2)^2+c^2*omega^2);
end
figure(1)
plot(OMEGA,C,'b',OMEGA,Ctheory,'ro-'); grid on;
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 results are not affected by the changes made to the initial conditions. The expression for C
does not depend on the initial conditions.

Exercise 3
% Exercise 3(a)
function LAB06ex2
clc
omega0 = 2; c = 0;
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) = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
end
figure(1)
plot(OMEGA,C,'ro-',OMEGA,Ctheory,'b'); grid on;

legend ('computed numerically','theoretical');

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

(a)

There is a sharp peak at the resonant frequency value at omega-naught-2, the solution skyrockets
until it equals 12.1
Max amplitude= 12.15; Frequency = 2; and 0 are both equal to 0.
The frequency yielding the maximal amplitude of the forced solution is 2. Compared to omega0,
they are both the same.
(b)
% Exercise 3 (b)
function LAB06ex1
clc
omega0 = 2; c = 0; 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-'); xlabel('t'); 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)]);

%---------------------------------------------------------------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 amplitude of the graph increases. This graph illustrates a vibration at resonance.

Exercise 4 :
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-'); xlabel('t'); 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)]);
%---------------------------------------------------------------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 ];

(a)
% Exercise 4(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);
C = 1/abs(omega0^2-omega^2);
A = 2*C*sin(0.5*(omega0-omega)*t);
figure(1)
plot(t,y,'b-',t,A,'r',t,-A,'g'); 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)]);
T=2*pi/omega
alfa=atan(c*omega/(omega0^2-omega^2))
hold on
yc=y-c*cos(omega*t-alfa);
plot(t,yc)
hold off
%---------------------------------------------------------------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 ];

(b)
computed amplitude of forced oscillation = 2.6137
theoretical amplitude = 1.3158
T=
3.4907
alfa =
0
Period (T) = 4 / |0-| = 4 / |2 - 1.8| = 62.83

(c)
Length of the beat: 2 / |0-| = 2 / |2 - 1.8| = 31.4159

(d)
% Exercise 4(d)
function LAB06ex1
clc
omega0 = 2; c = 0;
%omega = 1.9;
omega = 1.6;
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);
C = 1/abs(omega0^2-omega^2);
A = 2*C*sin(0.5*(omega0-omega)*t);
figure(1)
plot(t,y,'b-',t,A,'r',t,-A,'g'); 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)]);
T=2*pi/omega
alfa=atan(c*omega/(omega0^2-omega^2))
hold on
yc=y-c*cos(omega*t-alfa);
plot(t,yc)
hold off
%---------------------------------------------------------------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.9

Period (T): 4 / |0-| = 4 / |2-1.9| = 125.66

Length of the beat: 2 / |0-| = 2 / |2-1.9| = 62.83
The period of the fast oscillations decreased and the length of the beats increased.
2. = 1.6

Period (T) = 31.4159

Length of the beat = 15.70796

The period of rapid oscillations increased and the length of the beats decreased.
(e)
% Exercise 4(e)
function LAB06ex1
clc
omega0 = 2; c = 0; omega = 0.5;
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);
C = 1/abs(omega0^2-omega^2);
A = 2*C*sin(0.5*(omega0-omega)*t);
figure(1)
plot(t,y,'b-',t,A,'r',t,-A,'g'); 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)]);
T=2*pi/omega
alfa=atan(c*omega/(omega0^2-omega^2))
hold on
yc=y-c*cos(omega*t-alfa);
plot(t,yc)
hold off
%---------------------------------------------------------------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 ];

There are no beats present. This is so because there is no overlapping in the solution when omega
equals 0.5 or less. The initial condition and the current frequency are too far apart for
overlapping to occur.