EXERSIZE 1
clear ; all %this deletes all variables
omega0 = 4; c = 3; omega = 4.5;
param = [omega0,c,omega];
t0 = 0; y0 = 0; v0 = 0; Y0 = [y0;v0]; tf = 30;
options = odeset('AbsTol' 'relTol',1e-10, ,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 = 9; 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 ];
end
1
A)
Period of the forced oscillation is:
T=2π
ω=2π
4.5 ≈1.39
Numerical value of α:
α=π+tan−1
(
cω
ω0
2−ω2
)
, if ω ω
0<
π+tan−1
(
3(4.5)
4 4.5
2−2
)
