% damposcseries.m % % Function to plot x(t) for a damped oscillator % for a series of damping values % Takes initial conditions so that initial velocity is (nearly) zero % See also damposc.m % % RP 16 Oct. 2006; slight modifications Oct. 19, 2007 clear all close all t = 0.0:0.1:20; % array of time values A = 1.0; % "amplitude," in zero damping limit % gamma over 2 omega_0 -- make an array with logarithmically spaced % increments, using the "logspace" command. Could also do: % logg = -1.0:0.05:1.0; gw = 10.^logg; gw = logspace(-1,1,40); % Assign omega_0 = 1.0; N = length(gw); x = zeros(length(t),N); % initialize % Create the figure window figure; axis([min(t(:)) max(t(:)) -1 1]); hold on xlabel('time', 'FontWeight', 'bold'); ylabel('x', 'FontWeight', 'bold'); title('Damped Oscillator -- varying \gamma', 'FontWeight', 'bold'); text(10,-0.8, '\gamma / 2 \omega_0 = ', 'FontWeight', 'bold'); text(4,0.9, 'Initial amplitude A, initial velocity 0; \omega_0 = 1', ... 'FontWeight', 'bold'); for j=1:N, if (gw(j) < 1.0) % underdamped: g < 2*w0 w = sqrt(1.0 - gw(j)*gw(j)); % damped angular frequency alpha = atan(-1.0*gw(j)); Ampl = A / cos(alpha); x(:,j) = Ampl*exp(-gw(j)*t).*cos(w*t + alpha); % make previous plot gray if (j > 1) plot(t, x(:,j-1), 'color', [0.8 0.8 1.0], 'linewidth', 2.0); text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]); end plot(t, x(:,j), 'color', [0 0 0.5], 'linewidth', 2.0); text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold'); if (j==1) % pause, to explain... text(14, -0.4, 'paused...', 'FontWeight', 'bold') pause text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0]) end lastj = j; elseif (abs(gw(j)-1.0) < 0.001) % critically damped: g = 2*w0 if (j == lastj + 1) text(14, -0.4, 'paused...', 'FontWeight', 'bold') pause text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0]) text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]); end % for zero inital velocity, B = A*omega_0 x(:,j) = (A + A*t).*exp(-1.0*t); % make previous plot gray plot(t, x(:,lastj), 'color', [0.8 0.8 1.0], 'linewidth', 2.0); text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]); plot(t, x(:,j), 'k', 'linewidth', 3.0); text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold'); lastj = j; else % overdamped: g > 2*w0 if (j == lastj + 1) text(14, -0.4, 'paused...', 'FontWeight', 'bold') pause text(14, -0.4, 'paused...','FontWeight', 'bold','color', [1.0 1.0 1.0]) text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]); end beta = sqrt(gw(j)*gw(j) - 1.0); % for zero inital velocity, A1 = -A2*(gamma/2 - beta)/(gamma/2 + beta)' A2 = A / (1.0 - ((gw(j) - beta)/(gw(j) + beta))); A1 = -A2*(gw(j) - beta)/(gw(j) + beta); x(:,j) = A1*exp(-(gw(j) + beta)*t) + A2*exp(-(gw(j) - beta)*t); % make previous plot gray if (j > (lastj+1)) plot(t, x(:,j-1), 'color', [1.0 0.8 0.8], 'linewidth', 2.0); text(14,-0.8, num2str(gw(j-1)), 'FontWeight', 'bold','color', [1.0 1.0 1.0]); end plot(t, x(:,j), 'r', 'linewidth', 2.0); text(14,-0.8, num2str(gw(j)), 'FontWeight', 'bold'); end pause(0.5); end