matlab help

Exercise 1:

function LAB06ex1

clc

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

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

Exercise 2:

function LAB06ex2

omega0 = 3; c = 1;

OMEGA = 2:0.02:4;

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) = ??; % FILL-IN

end

figure(2)

plot(??); 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 ];