% @(#) 6.728 Homework Problem 3.1 -- Matlab Code % ps3p1.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 % .m file written by Mark Schweizer % modified by V. Anant % 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 = (L^2*2*pi)^(-1/4) * exp(-(x / (2*L)).^2); 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; 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; % 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); title('P3.1: Standard Deviation in x-space'); xlabel('Time [fs]'); ylabel('\Delta x [nm]'); 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;