1 / 11100%
Mat 275 MATLAB LAB #6
Name: Brandon Grebe
Lab Day and Time: Wednesday
10:30
Instructor: Chowell Puente
Exercise #1
Part A
By running the given function file we can see the graph of the solution as
well as find the period of this forced oscillation.
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)]);
%----------------------------------------------------------------
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 ];
Looking at this graph we can see that the period is going to be somewhere
between 4 and 5.
We can find the exact period by using the following equation and knowing
omega = 1.4:
T = 2
2𝜋
𝜔
=
2𝜋
1.4
= 𝟒. 𝟒𝟗
To find the numerical angle 𝛼 we can use the arctan equation. Knowing c =
1, 𝜔 = 1.4, and 𝜔0 = 2.
𝛼 = arctan(c
𝜔
𝜔0
2
−𝜔
2
) = arctan (
1
(
1.4
)
2
2
−1.4
2
) = 0.601 radians
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));
yc=y-Ctheory*cos(omega*t-alpha);
figure(2)
plot(t,yc,'b-');
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 ];
Looking at this graph we can confirm that this solution is exponentially
decreasing, we can see that it dies down around t = 10 which agrees with the
particular solution where at t = 10 the solution continues to oscillate at
the same level.
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(k) = 1/sqrt((omega0^2-
omega^2)^2+(c*omega)^2); % FILL-IN
end
figure(2)
plot(OMEGA,C,OMEGA,Ctheory,'ro'); 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 ];
Part B
To determine the maximum value of C we can simply zoom in on the graph to
find when C is greatest. We can see this occurs around 𝜔 = 1.86 and C = 0.5162
Part C
Part 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)]);
alpha=atan(c*omega/(omega0^2-omega^2));
yc=y-Ctheory*cos(omega*t-alpha);
figure(2)
plot(t,yc,'b-');
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 ];
Running the solution with the omega value we found in the last problem we are
presented with this graph. We can see that the oscillations have indeed
increased as we expected. Zooming in we can see that the amplitude is 0.5. If
omega is increased we expect to see a larger amplitude, the same goes for a
smaller omega, the amplitude will be smaller.
Part E
Exercise #3
Part A
The maximal amplitude approaches infinity as omega approaches 2.
However, we can see that the computed maximum amplitude is roughly 12. The
value of 2 for omega would yield the greatest amplitude.
Part 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-'); 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));
yc=y-Ctheory*cos(omega*t-alpha);
figure(2)
plot(t,yc,'b-');
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 ];
This graph confirms what we saw in the last exercise. We can see that the
oscillations are increasing linearly without stopping, these oscillations
would approach infinity of course if we had a larger window.
Exercise #4
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;
C=1/abs(omega0^2-omega^2);
A=2*C*sin((omega0-omega)*t/2);
hold on
plot(t,A,'r-',t,-A,'g-');
%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));
%yc=y-Ctheory*cos(omega*t-alpha);
%figure(2)
%plot(t,yc,'b-');
%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 ];
Part B
The period of fast oscillation can be found analytically by using the
following equation:
𝑇 =
4𝜋
𝜔0 + 𝜔
=
4𝜋
2 + 1.8
= 𝟑. 𝟑𝟎𝟕
Part C
The length of each beat can be found by using the following equation:
4𝜋
𝜔0 𝜔
=
4𝜋
0.2
= 20𝜋
Thus the length of each beat is 10 𝜋=31.42
Part D
Omega = 1.9:
Fast Oscillation:
3.22
Length of Beat:
62.83
Omega = 1.6:
Fast Oscillation:
3.49
Length of Beats:
15.707
Part E