%%%%%% Problem Set 2 %%%%%%%%%%% % Problem 2.2 % V. Anant (modification of original code by B. Goldberg) % further modified W.M. Kaminsky 9/24/2006 % - change units to nanometers throughout % - make prettier graph (annotation, xlim, higher res analytic % expression) % - fixed mislabeling of standard deviation as variance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % See the notes on the 6.728 website for information % on taking a discrete fourier transfom % preamble clear close all format compact disp('===================='); disp('Problem Set 2 output'); disp('===================='); % Define an x vector x = [-20:0.02:20].'; % Define an L vector L = [0.5 1 2].'; % Define a matrix of Psi_L(x) Psi(:,1) = (pi*L(1)^2)^(-1/4).*exp(-x.^2/(2*L(1)^2)); Psi(:,2) = (pi*L(2)^2)^(-1/4).*exp(-x.^2/(2*L(2)^2)); Psi(:,3) = (pi*L(3)^2)^(-1/4).*exp(-x.^2/(2*L(3)^2)); %%%% Part b %%%%% deltax = 0.02; N = (20 - -20)/deltax + 1; deltaq = 2*pi/(N*deltax); qmax = pi*(N-1)/(deltax*N); qmin = -qmax; % Create q vector q = [qmin:deltaq:qmax].'; % Create phi matrix phi = exp(i*(x*q.')); % Calculate the fourier transform A(q) A = phi'*Psi*deltax; clf; plot(q, abs(A(:,1)).^2,'x',q, abs(A(:,2)).^2,'o',q, abs(A(:,3)).^2,'+'); %%%%%%%% Part c %%%%%%%%% disp('Part (c): A_L normalization'); disp('-----------------------'); fprintf('Sum should be 1 for all L''s\n\n'); Sum = [ sum(A(:,1)'*A(:,1)*deltaq/(2*pi)); sum(A(:,2)'*A(:,2)*deltaq/(2*pi)); sum(A(:,3)'*A(:,3)*deltaq/(2*pi)); ] % Compare against the analytic expression for A(q) qa = [-5:0.002:5].'; %% Essentially the entire nonzero part of the Fourier %% transform happens between q = -5/nm to 5/nm %% To get a nice Gaussian shape, we plot it out at %% 10x higher resolution than in q A_analytic(:,1) = (4*pi*L(1)^2)^(1/4)*exp((-qa.^2*L(1)^2)/2); A_analytic(:,2) = (4*pi*L(2)^2)^(1/4)*exp((-qa.^2*L(2)^2)/2); A_analytic(:,3) = (4*pi*L(3)^2)^(1/4)*exp((-qa.^2*L(3)^2)/2); hold on; plot(qa, abs(A_analytic).^2); title('Fourier Transform of \psi(x) = [1/(\pi L^2)]^{1/4} e^{-x^2 /(2L^2)}'); xlim([-5 5]) %% Again, essentially the entire nonzero part of the Fourier transform happens between q = -5/nm to 5/nm xlabel('q [nm^{-1}]'); ylabel('|A(q)|^2'); legend('L = 0.5nm', 'L = 1nm', 'L = 2nm'); annotation1 = annotation('textbox',... 'Position',[0.6175 0.7014 0.245 0.06968],... 'FitHeightToText','off',... 'String',{'Solid lines denote analytical expression (4\piL^2)^{1/4} e^{-q^2 L^2 /2}','','x,o, and + signs denote numerical discrete Fourier transform.'}); hold off; %%%% Part d %%%%% disp('Part (d): variance of x,q'); disp('-----------------------'); mean_x = [ sum(Psi(:,1)'*(x.*Psi(:,1))); sum(Psi(:,2)'*(x.*Psi(:,2))); sum(Psi(:,3)'*(x.*Psi(:,3)))]*deltax; meansq_x = [ sum(Psi(:,1)'*(x.^2.*Psi(:,1))); sum(Psi(:,2)'*(x.^2.*Psi(:,2))); sum(Psi(:,3)'*(x.^2.*Psi(:,3)))]*deltax; std_x = sqrt(meansq_x - mean_x.^2); mean_q = [ sum(A(:,1)'*(q.*A(:,1))); sum(A(:,2)'*(q.*A(:,2))); sum(A(:,3)'*(q.*A(:,3)))]*deltaq/(2*pi); meansq_q = [ sum(A(:,1)'*(q.^2.*A(:,1))); sum(A(:,2)'*(q.^2.*A(:,2))); sum(A(:,3)'*(q.^2.*A(:,3)))]*deltaq/(2*pi); std_q = sqrt(meansq_q - mean_q.^2); uncertainty = abs(std_q.*std_x); disp('Standard Deviation (Uncertainty) of x') disp(std_x); disp('Standard Deviation (Uncertainty) of q') disp(std_q); disp('Uncertainty Product') disp(uncertainty)