% fourierdemo.m % % Program to illustrate Fourier Series decomposition of a vibrating string. % Can also show oscillations (normal mode decomposition) as a function of % time. % % Raghuveer Parthasarathy % Nov. 2007 % % Options (functions) % 1. Triangle % y(x) = kx, x < L/2 % kL - kx, x >= L/2 % Fourier: sine terms Bn = 4 *sin(j*pi/2)/ (n pi)^2 (odd n) % Bn = 0 (even n) % 2. Weird rectangle function % y(x) = 0 for x < L/4 % = 1 for L/4 <= x < 3L/5 % = -1 for 3L/5 <= x < 7L/8 % = 0 for 7L/8 <= x % 3. Gaussian % y(x) = exp(-(x-3L/4).^2 / 2 / 0.05 / 0.05); % % Raghuveer Parthasarathy, 5 Nov. 2006 clear all close all % ---------------- disp('* Fourierdemo.m -- RP Nov. 2007') disp(' '); disp('Functions: ') disp(' (1) Triangle'); disp(' (2) Weird rectangle function'); disp(' (3) Gaussian'); opt = input('Enter option: '); disp(' '); N = input('Number of terms (max. n value) of Fourier Series: '); % number of terms in series oscopt = logical(input('Display oscillations? (1==yes): ')); % ------------------------------------------ x = 0.0:0.005:1.0; % x, in units of L Nx = length(x); % Define the functions switch opt case 1 %triangle y = x.*(x <= 0.5) + (1.0-x).*(x > 0.5); % y, in units of k case 2 % weird rectangle y = zeros(1, Nx); for j=1:Nx, if (x(j) < 0.25), y(j) = 0.0; elseif (x(j) < 0.6), y(j) = 1.0; elseif (x(j) < 0.875), y(j) = -1.0; else y(j) = 0.0; end end case 3 % Gaussian y = exp(-(x-3/4).^2 / 2 / 0.05 / 0.05); otherwise disp('Error! Invalid option for funtion. Hit Control-C'); pause end % Fourier components (cosine amplitudes are zero) B = zeros(1,N); switch opt case 1 % triangle % Analytically determine Fourier components for j=1:N, B(j) = 4 *sin(j*pi/2)/ j / j/ pi / pi; ycomp(j,:) = B(j)*sin(j*pi*x); end case {2 , 3} % weird rectangle, or Gaussian % computationally determine Fourier components -- I'm being lazy! for j=1:N B(j) = 2 * sum(y.*sin(j*pi*x)*mean(diff(x))); ycomp(j,:) = B(j)*sin(j*pi*x); end end % partial sum after N terms Ntermy = zeros(size(y)); for j=1:N, Ntermy = Ntermy + ycomp(j,:); end % plot Fourier sum after N terms figure(1) plot(x, y, 'k-', 'LineWidth', 2.0) xlabel('x / L', 'Fontweight', 'bold'); ylabel('y / k', 'Fontweight', 'bold'); title('y(x), and N term Fourier Sum', 'Fontweight', 'bold'); hold on plot(x, Ntermy, 'b-', 'LineWidth', 2.0) legend('y(x)', 'N term Fourier series', 'Location', 'Northwest'); legend('boxoff') axis auto % Plot each component figure(2) plot(x, y, 'k-', 'LineWidth', 2.0) xlabel('x / L', 'Fontweight', 'bold'); ylabel('y / k', 'Fontweight', 'bold'); title('y(x), and each Fourier component', 'Fontweight', 'bold'); hold on for j=1:N, plot(x, ycomp(j,:), 'color', [0 1 0]*j/N/4 + [0 0.75 0], 'LineWidth', 2.0) end axis auto a = axis; axis([a(1) a(2) (a(3)-0.1) (a(4)+0.1)]) % Plot Fourier amplitudes figure(3) plot(B, 'ko-', 'LineWidth', 2.0) xlabel('n', 'Fontweight', 'bold'); ylabel('B_n', 'Fontweight', 'bold'); title('Fourier amplitudes B_n', 'Fontweight', 'bold'); grid on if oscopt==1, % Illustrate oscillations disp(' '); disp('Illustration of oscillations.'); disp('Press enter to continue.'); pause disp(' ** MAKE FIGURE WINDOW BIGGER; Then Press Enter.'); figure; pause w0 = 2*pi*0.2; % 0.2 Hz "natural frequency" % Frequencies of the normal modes omega = zeros(1,N); for j=1:N, omega(j) = j*w0; end % calculate components and sum at each time interval k = 1; dt = 0.1; for t=0:dt:6*pi/omega(1), yfourier(k,:) = zeros(size(y)); cla for j=1:N ycomptime(k,j,:) = ycomp(j,:)*cos(omega(j)*t); yfourier(k,:) = yfourier(k,:) + ycomp(j,:)*cos(omega(j)*t); end k = k+1; end Nt = k-1; % Illustration of oscillations axis([-0.05 1.05 min(min(yfourier)) max(max(yfourier))]); hold on grid on title('Oscillations of non-normal modes.', 'Fontweight', 'bold'); xlabel('x', 'Fontweight', 'bold'); ylabel('y (Fourier sum)', 'Fontweight', 'bold'); dt = 0.05; for k = 1:Nt, % clear previous plot cla for j=1:N plot(x, squeeze(ycomptime(k,j,:)), 'color', ... [1 1 1]*0.7*(1.0 - (j/N)), 'Linewidth', 2.0); end plot(x, yfourier(k,:), 'r', 'linewidth', 2.0); if (k==1) disp('Press Enter to begin...'); pause end pause(0.1); end end