% @(#) 6.728 Homework Problem 3.3 -- Matlab Code % ps3p3_movie.m file % @(#) This program computes decomposes an the initial state of a particle % in a box state into the M lowest energy eigenstates (i.e., cosines and % sines). It then makes a movie of the time evolution. % @(#) See Also: Eigenfunction Expansion Tutorial % preamble 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 = 1; % 1D Box [-1/2, 1/2] is of unit length n_max = 20; % Create our x vector dx= 0.5E-3; x = [-0.5:dx:0.5]'; % Create the column vectors psi. psi = zeros(2001,1); psi(876:1126) = sqrt(8); % Initial state: psi(x,0) = \sqrt{\frac{8}{L}} if |x| \leq L/16 % else 0. x(876) = -0.5+875*0.0005 = -0.0625 = -L/16 as L = 1 % x(1126) = -0.5+1125*0.0005 = 0.0625 = L/16 as L = 1 psi_pdf = conj(psi).*psi; % Create matrix of basis functions phi where energy level n is incremented % column by column and location x is incremented row by row. % % phi_1(x(1)) phi_2(x(1)) ... phi_{n_max}(x(1)) % phi = phi_1(x(2)) phi_2(x(2)) ... phi_{n_max}(x(2)) % | | \ | % phi_1(x(2001)) phi_2(x(2001)) ... phi_{n_max}(x(2001)) % % in our case % phi_{n odd} = sqrt{2/L) * cos(n*pi*x/L) % phi_{n even} = sqrt(2/L) * sin(n*pi*x/L) % phi = zeros(2001,n_max); % allocate space for phi matrix N = ones(2001,1)*[1:n_max]; X = x*ones(1,n_max); if rem(n_max,2) == 0 % that is, if M is even phi(:,1:2:n_max-1) = sqrt(2/L)*cos(N(:,1:2:n_max-1).*pi.*X(:,1:2:n_max-1)/L); % odd columns are cosines phi(:,2:2:n_max) = sqrt(2/L)*sin(N(:,2:2:n_max).*pi.*X(:,2:2:n_max)/L); % even columns are sines else % that is, if M is odd phi(:,1:2:n_max) = sqrt(2/L)*cos(N(:,1:2:n_max).*pi.*X(:,1:2:n_max)/L); % Again, odd columns are cos phi(:,2:2:n_max-1) = sqrt(2/L)*sin(N(:,2:2:n_max-1).*pi.*X(:,2:2:n_max-1).L); % even columns are sines end a = dx * phi' * psi; % Decomposition of initial state psi into the % M lowest energy levels (cosines and sines) a_pdf = abs(a).^2; % Create our energy and time matrices m = Me; % we'll say the particle is an electron Eg = (hbar*pi)^2 / (2*m*L^2); % Ground state energy (n = 1) varrho = m*L^2/h; % We'll measure time as fractions of the t = varrho*[0:0.01:1]'; % reconstruction time mL^2 / h Ecol = [(hbar*pi)^2/(2*m*L^2)]*[1:n_max].^2' ; E = Ecol*ones(1,length(t)); T = ones(n_max,1)*t'; eta = exp(-i*E.*T/hbar); a_t = (a*ones(1,length(t))).*eta; psi_t = phi*a_t; psi_t_pdf = conj(psi_t).*psi_t; % Make movie for j = 1:length(t); plot(x,psi_t_pdf(:,j)) text(0.2,8,['Time t = ',num2str(t(j)/varrho),'mL^2/h']) xlabel('Position x (in units of box length L)') ylabel('|\Psi(x,t)|^2 (probability per unit length)') axis([-0.5 0.5 0 max(max(psi_t_pdf))]); M(j) = getframe; end