% Ncoupledosc.m % % Program to illustrate normal modes of N coupled oscillators % Will also illustrate "redundancy" of solutions for "modes" n>N % Loaded string % % Raghuveer Parthasarathy % October 28, 2006 clear all close all N = round(input('Enter the number of oscillators (4): ')); % Will consider 2*N+1 modes, though of course "unique" modes are only 1:N w0 = 2*pi*0.2; % 0.2 Hz "natural frequency" C = 1; % same amplitude for all modes % beads on string, separated by l=1; string goes from x = 0 to N+1 % Frequencies of the modes omega = zeros(1,2*N+1); for j=1:2*N+1, omega(j) = 2*w0*sin(j*pi / 2 / (N+1)); end % Amplitude array; dim 1 = p, dim 2 = n % Also include endpoints (zero amplitude), so point p is at index p+1 A = zeros(N+2,2*N+1); for n = 1:2*N+1, % modes for p=1:N+2, % particles A(p,n) = C*sin((p-1)*n*pi/(N+1)); end end % Will also plot continuous sine function along x: x = 0:0.05:N+1; Nx = length(x); y_cont = zeros(2*N+1,Nx); for n=1:2*N+1, y_cont(n,:) = C*sin(x*n*pi/(N+1)); end % ------------------------ % Options disp('coupledosc.m -- R. Parthasarathy 28Oct06'); disp('** Options **'); disp(' (1) Illustration of normal modes at t=0.'); disp(' (2) Illustration of oscillations of normal modes.'); disp(' (3) Illustration of oscillations of normal modes for "modes" n>N.'); opt = input('Enter the option number: '); switch opt case 1 disp('Illustration of normal modes at t=0.'); % Illustration of normal modes at t=0 % Initial conditions: \delta=0 (zero initial velocity) figure; subplot(N+1,1,1); for j=1:N+1, subplot(N+1,1,j); plot(x, y_cont(j,:), 'g-', 'Linewidth', 2.0); hold on plot([0:N+1], A(:,j), 'bo', 'LineWidth', 2.0); plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black grid on axis([-0.1 N+1+0.1 -1.2 1.2]); if (j==1) title('normal modes at t=0', 'Fontweight', 'bold'); end if (j==N+1) xlabel('x', 'Fontweight', 'bold'); end ylabel('y', 'Fontweight', 'bold'); end case 2 disp('Illustration of oscillations of normal modes.'); disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.'); figure; pause % Illustration of oscillations subplot(N,1,1); dt = 0.2; for t=0:dt:4*pi/omega(1), for j=1:N, subplot(N,1,j); % clear previous plot cla plot(x, y_cont(j,:)*cos(omega(j)*t), 'color', [0.7 0.7 0.7], 'Linewidth', 2.0); hold on plot([0:N+1], A(:,j)*cos(omega(j)*t), 'bo-', 'LineWidth', 2.0); plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black grid on axis([-0.1 N+1+0.1 -1.2 1.2]); if (j==1) title('Oscillations of normal modes.', 'Fontweight', 'bold'); end if (j==N+1) xlabel('x', 'Fontweight', 'bold'); end ylabel('y', 'Fontweight', 'bold'); end if (tN.'); disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.'); figure; pause % Illustration of oscillations -- aliased mode number (n > N) subplot(N+1,2,1); dt = 0.2; for t=0:dt:4*pi/omega(1), % for j=1:(2*N+1), for j=1:(2*N+1), ycont_t(j,:) = y_cont(j,:)*cos(omega(j)*t); % continuous function ypart_t(j,:) = A(:,j)*cos(omega(j)*t); % discrete particle positions end for j=1:(2*N+1), plotno = (j<(N+2))*(2*j-1) + (j>(N+1))*((2-2*N)*j/(N-1) + 4*(N+1)); % which subplot element subplot(N+1,2,plotno); % clear previous plot cla plot(x, ycont_t(j,:), 'color', [0.7 0.7 0.7], 'Linewidth', 2.0); hold on plot([0:N+1], ypart_t(j,:), 'ro-', 'LineWidth', 2.0); plot([0 N+1], [0 0], 'ko', 'LineWidth', 2.0); %ends in black grid on axis([-0.1 N+1+0.1 -1.2 1.2]); if (j==N+1) xlabel('x', 'Fontweight', 'bold'); end ylabel('y', 'Fontweight', 'bold'); fs = sprintf('Mode n=%d', j); title(fs, 'Fontweight', 'bold'); end if (t