% @(#) 6.728 Homework Problem 3.1 -- Matlab Code % ps3p1mod.m file % @(#) This program computes a DFT by hand. It then performs a time evolution. % @(#) See Also: Eigenfunction Expansion Tutorial % 6.728 -- Problem 3.1, part a % .m file originally written by Mark Schweizer % modified by V. Anant % Further modified 9/26/2006 by W.M. Kaminsky % Corrected psi to agree with textbook/prob set convention (line 38 ff.) % Made new code for analytical solution as old code seemed to disagree % with numerical solution for all times t > 0. (lines 90-107) % preamble diary ps3p1.out % keep a log of our output diary on; clear all; % Clear workspace clc; % Clear the command window close all; % Get rid of old figures format compact; % suppresses extra lines % Defines some standard constants used in quantum mechanics and % statistical mechanics. h = 6.62607E-34; % Planck constant [J*s] hbar = h/(2*pi); % Rationalized Planck constant [J*s] Me = 9.1094E-31; % Electron rest mass [kg] Mp = 1.672623E-27; % Proton rest mass [kg] Mn = 1.674929E-27; % Neutron rest mass [kg] c = 2.99792458E8; % Speed of light in vacuum [m/s] L = .5 *1E-9; % Gaussian Width Paramters % Create our x vector dx=.2E-9; x = [-20E-9:dx:20E-9]'; % Create the column vectors psi. psi = (pi*L^2)^(-1/4) * exp(-[x.^2]./[2*L.^2]); %%%% WMK 9/26/2006 %%%% % Code previously had psi = (L^2*2*pi)^(-1/4) * exp(-(x / (2*L)).^2); % which implies abs(psi).^2 = (L^2*2*pi)^(-1/2) * exp(-[x.^2]./[2*L.^2]); % This is a correctly normalized Gaussian with mean 0 and stand dev L % However, PS3 lists psi as % psi = (pi*L^2)^(-1/4) * exp(-[x.^2]./[2*L.^2]); % --> abs(psi).^2 = (pi*L^2)^(-1/2) * exp(-[x.^2]./[L.^2]); % which is a correctly normalized Gaussian with mean 0 and stand dev % L/sqrt(2). Personally, I think the former convention makes more sense % than the latter since abs(psi).^2 is the actual probability % distribution. Nevertheless, I corrected m-file to agree with textbook % and problem set convention. %%%%%%%%%%%%%%%%%%%%%%% psi_pdf = conj(psi).*psi; % Here we set up the q vector. "dq" is determined by the periodicity % of the DFT, as is the Max/Min q available. N=size(x,1); % N = total number of points in x dq=2*pi/(N*dx); q=dq*[-(N-1)/2:(N-1)/2]'; % Compute our basis matrix and determine the expansion coefficients phi = exp(i*(x*q')); a = dx * phi' * psi; a_pdf = abs(a).^2; % Create our energy and time matrices % Note, that in all that follows, I've again combined a series of % vectors into one matrix. For example, the energy matrix defined % below consists of columns of energy vectors as defined in the % tutorial. Each column corresponds to a different time instance. tau = Me*L^2/hbar t = tau*[0 2 5 10 15]'; E = hbar^2/(2*Me) * q.^2; eta = exp(-i* E/hbar*t.'); a_t = (a*ones(size(t'))).*eta; psi_t = phi *a_t * dq /(2*pi); psi_t_pdf = conj(psi_t).*psi_t; % Plot |psi|^2 for each t. figure; colormap('white'); plot(x/1E-9,psi_t_pdf); % make sure cyclelines.m is in the same directory as this file cyclelines; axis([-10 10 0 max(max(psi_t_pdf))]); % Create a vector with the analytical solution m = Me; %%%%% WMK 9/26/2006 %%%%%%%% %% Make my own code for plotting analytical solution since it seems to me %% that the original code, namely % %psia = ones(size(x))*(((2*pi*L^2*(1+i*t*hbar/(L^2*2*m)).^2).^(-1/4)).') ... % .*exp(-(1/(4*L^2))* (x.^2)*(1./(1+i*t*hbar/(2*m*L^2))).'); %psia_pdf = conj(psia).*psia; % %% has some sort of error as it only agrees well with numerical solution at %% time = 0. X = x*ones(1,length(t)); % Make matrices (X,T) with x incremented row by T = ones(length(x),1)*t'; % row and t incremented column by column var_t = (0.5*L^2)*[1 + (hbar/(m*L^2))^2 * T.^2]; %% Textbook Equ 4.20 psia_pdf = [(2*pi*var_t).^(-1/2)].*exp(-[X.^2]./(2*var_t)); %% Equ. 4.21 %%%% END WMK 9/26/2006 %%%%%%%% % Plot the analytical solution on the same graph hold on; plot(x/1E-9,psia_pdf,'x'); legend('t=0','t=2 tau','t=5 tau','t=10 tau','t=15 tau','Analytic solutions'); title('P3.1: Gaussian Wave Functions in Time'); xlabel('X [nm]'); ylabel('\Psi* \Psi [m^{-1}]'); prob_sums_x= diag(psi_t'*psi_t)*dx mean_x = (diag(psi_t'*diag(x,0)*psi_t)*dx) mean_sq_x = (diag(psi_t'*diag(x.^2,0)*psi_t)*dx) std_dev_x=sqrt(mean_sq_x - mean_x.^2) prob_sums_q= diag(a_t'*a_t)*dq/(2*pi) mean_q = (diag(a_t' * diag(q,0) * a_t ) * dq/(2*pi)) mean_sq_q = (diag(a_t' * diag(q.^2,0) * a_t )* dq/(2*pi)) std_dev_q=sqrt(mean_sq_q - mean_q.^2) hup = std_dev_x.*std_dev_q figure; % we take the real part because the calculations above % generates a small imaginary part due to machine precision % errors. plot(t/1e-15,real(std_dev_x)/1E-9,'-',t/1e-15,1E9*sqrt(var_t(1,:)),'x'); title('P3.1: Standard Deviation in x-space'); xlabel('Time [fs]'); ylabel('\Delta x [nm]'); legend('Numerical solution','Analytical solution') figure plot(t/1e-15,real(std_dev_q)*1E-9); title('P3.1: Standard Deviation in q-space'); xlabel('Time [fs]'); ylabel('\Delta q [nm^{-1}]'); diary off;