MAT 275 - MATLAB #6_v2
% Solution 1
On analyzing the behavior of the forced oscillations using the figure L6a, we can see that it is
identical to a sinusoidal waveform. Therefore, we can calculate the time-period using the
method of visual inspection. The time-period comes out to be approximately equal to 2.5 secs.
Also on using the equation L6.4a α = 0.635 radians.
% (B)
The updated function is attached below -
function LAB06ex1
omega0 = 3; c = 1; omega = 2.4;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 50;
options = odeset('AbsTol',1e-10,'RelTol',1e-10);
% Equation to calculate “alpha”
alpha = atan(omega/(omega0^2 - omega^2));
[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;
% Complimentary solution is defined
Yc = y - C.*cos(omega.*t - alpha);
figure(2)
plot(t,Yc,'r-'); ylabel('Yc'); grid on;
disp(['computed amplitude of forced oscillation = ' num2str(C)]);
Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
disp(['theoretical amplitude = ' num2str(Ctheory)]);
%--------------------------------------------------- %
end
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 ];
end
Complimentary Solution
From the above figure, we can see that the resulting figure is like exponentially decaying
oscillations.
The sinusoidal terms (sine and cosine) gives the oscillations but the exponential term is the
reason for decaying amplitude.
% Solution 2
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;
% New Ctheory calculated here
Ctheory(k) = 1/sqrt((((omega0^2)-(omega^2))^2)+(c*omega)^2); % FILL-IN
end
figure(2)
% Updated plot function
plot(OMEGA,C, 'k-', 'LineWidth', 2, OMEGA,C, 'ro'); grid on; % FILL-IN
xlabel('\omega'); ylabel('C');
end
%---------------------------------------------------------
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 ];
end
% (B)
Fig - C vs ω
Using the above figure, we can calculate the approximate value of Omega for which C is
maximum is:
ω = 2.92
C = 0.338
% (C)
To calculate the value of ω for which C is maximum, we can use the following equation -
𝐶(𝑤) = 1
((𝜔𝑜
2− 𝜔2)2+ 𝑐2𝜔2)1/2
At maxima, the differentiation of C(ω) w.r.t ω should be equal to zero. Using this fact :
−4𝜔(𝜔𝑜
2− 𝜔2) + 2𝑐2𝜔 = 0
𝜔 = (𝜔𝑜
2−𝑐2
2)1/2
Putting ωo as 3 and c = 1, we can calculate ω -
ω = 2.91
% (D)
The above value of ω is obtained for the case when “c” is maximum. Therefore, at this value the
amplitude of the forced oscillations would also be maximum. Any other value of ω would result
in the reduction of amplitude.
% (E)
Even after modifying the initial conditions {y(0) = 1 and y’(0) = 1}, the curve doesn’t change at
all -
Fig – C vs Omega
% Solution 3
% (A)
For resonance to occur the value of c should be equal to 0.
Using the above value, we can get the maximum value of the amplitude.
The maximum value of amplitude is 8, this occurs when ω is equal to 3:
The value for ω and ωo comes out to be same {case of resonance}
% (B)
Curve obtained is shown below -
Fig – Decreasing Sinusoidal
The curve is similar to exponentially decreasing sinusoidal wave.
% Solution 4
function LAB06ex1
omega0 = 2; c = 0; omega = 1.8;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 100;
alpha = atan(omega/(omega0^2 - omega^2));
options = odeset('AbsTol',1e-10,'RelTol',1e-10);
[t,Y] = ode45(@f,[t0,tf],Y0,options,param);
y = Y(:,1); v = Y(:,2);
% Calculating correct “C” here -
if omega > omega0
C = 1/(omega^2 - omega0^2);
else
C = 1/(omega0^2 - omega^2);
end
% The envelope function A is defined below
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;
%yc = y - C.*cos(omega.*t - alpha); % Don’t need the complimentary solns
%figure(2)
%plot(t,yc,'r-'); ylabel('yc'); grid on;
disp(['computed amplitude of forced oscillation = ' num2str(C)]);
Ctheory = 1/sqrt((omega0^2-omega^2)^2+(c*omega)^2);
disp(['theoretical amplitude = ' num2str(Ctheory)]);
end
%------------------------------------------------------------
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 ];
end
% (A)
Fig – Envelope
% (B)
From the document the period of the fast oscillations is equal to:
T = 4π/ (ω + ωo) ωo = 3 & ω = 2.8 rad/s
T = 2.32sec
The time-period obtained from the graphs as well as shown below:
Fig – Time Period (1)
Fig – Time Period (2)
% (C)
As per the document the length of the beat is equal to:
L = 2π/ (ωo - ω) ωo = 3 & ω = 2.8 rad/s
L = 31.42
From the graph a similar value is obtained as shown below:
Fig – Beat Length
% (D)
Changing value of Omega to 2.9:
Fig – Omega = 1.9
Time period of the beats: 2.12s Length of the beats: 62.83
Changing the value of omega to 2.6 rad/s
Fig 10 – Omega = 1.6
Time period of the beats: 2.24s
Length of the beats: 15.7
From the two figures above we can infer that as the value of omega gets closer to the value of
Omega0, the time period decreases when compared to the previous value of omega (2.8 rad/s).
Though the time period decreases but we can see that there is a considerable increase in the
length of the beats as the value of omega approaches Omega0.
Further, on decreasing the value of omega and moving it away from the value of Omega0, the
time period increases and the length of the beats decreases by a considerable amount.
% (E)
At Omega = 1.5 the phenomenon of the beats will not be observed.
The phenomenon of having a very rapidly varying function combined with a slowly varying
function leads to beats. Such scenario is created when the value of Omega is close enough to
the value of Omega0. The given value of Omega (=1.5) is quite far from Omega0. Thus beats will
not be observed for this value of Omega.
