% % Particle scattering from a potential well. We time evolve a % starting wave packet. The calculations are done in a simple % brute force way: we set up a large periodic box of size L, % discretize it, setup the kinetic and potential operators, get H, % diagonalize it, and then project the intitial state onto % eigenstates and time evolve directly in the eigenbasis. So it is % exact but expensive. % % Just set key constants to one for now although we do set constants % properly below. clear hbar=1; m=1; % Size of box and then we discretize the x axis from -L/2 to L/2 with nx % points. dx is grid spacing. L=200; nx=500; x=linspace(-L/2,L/2,nx)'; dx=x(2)-x(1); fprintf('L=%g nx=%d dx=%g\n',L,nx,dx); % Potential barrier/well : Maximum height/depth is V0 and width a % centered at origin so from -a/2 to a/2 a = 2; fprintf('Barrier width a=%g from -a/2 to a/2 [-%g to %g]\n',a,a/2,a/2); V0 = 1; V = V0*exp(-(x/(a/2)).^10); fprintf('Maximum potential height = %g\n',V0); % Set up Hamiltonian with PBC % direct Fourier representation of Laplacian (D2) kvec = 2*pi/(nx*dx)*[-nx/2:nx/2-1]'; cD2 = ifft(fftshift(-kvec.^2)); for j=1:nx D2(:,j) = circshift(cD2,j-1); end D2 = (D2+D2')/2; T = -hbar^2/(2*m)*D2; H = T+diag(V); [psin,E] = eig(H); E = diag(E); psin = psin/sqrt(dx); fprintf('First 5 eigenvalues: %+g\n',E(1:5)); % Initial state psi0=psi(x,t=0). % A Gaussian packet with width sig and momentum hbar*k0: sig = 20; k0 = sqrt(2*m*.5)/hbar; psi0 = exp(-(x+L/4).^2/sig^2).*exp(i*k0*x); % Normalize psi0 and then compute the expansion coefficients c(:) in terms % of the eigenfunctions. psi0 = psi0/sqrt(sum(abs(psi0).^2)*dx); c = psi0'*psin*dx; c = c'; E0 = sum(abs(c).^2 .* E); fprintf('Initial packet has energy E0 = %g with sigE = %g\n',... E0,sqrt(sum(abs(c).^2 .* E.^2)-E0^2)); if (V0>0) idxV0=find(abs(E-V0)==min(abs(E-V0))); idxmax=find(abs(c)==max(abs(c))); fprintf('Prob of energy E=%g is %g\n',E(idxmax),abs(c(idxmax))^2); fprintf('Prob of energy E=%g is %g\n',E(idxV0),abs(c(idxV0))^2); end % Get some maximum values to set scale for plots maxp = max(abs(psi0).^2)*1.5; maxr = sqrt(maxp); % We will to nt time steps from t=0 to t=tmax % Set to max to be the time for the center of wave packet moving % with group velocity hbar*k0/m to go a distance of L/2 (half box) fs=20; counter=0; nt=90; tmax=(L/2)/(hbar*k0/m); %30*m*sig^2/(hbar*2*pi); % x ranage of plots xmin=-L/2; xmax=L/2; for t=0:tmax/nt:tmax counter=counter+1; % compute psit = psi(x,t) using time-evolved eigenbasis expression psit = psin*(c.*exp(-i*t*E/hbar)); % make 3 stacked plots and write to file clf subplot(3,1,1) plot(x,real(psit),'b'); set(gca,'fontsize',fs,'fontname','Arial') ylabel('Real \Psi') grid axis([xmin xmax -maxr maxr]); title([sprintf('t = %.3f ',0.5*t/tmax) 'L*m/(hbar*k0)']) subplot(3,1,2) plot(x,imag(psit),'r'); set(gca,'fontsize',fs) ylabel('Imag \Psi') grid axis([xmin xmax -maxr maxr]); subplot(3,1,3) plot(x,abs(psit).^2,'k'); set(gca,'fontsize',fs,'fontname','Arial') xlabel('x','fontname','Arial'); ylabel('|\Psi|^2','fontname','Arial') grid axis([xmin xmax 0 maxp]); fname = sprintf('frame_%04d.jpg',counter-1); fprintf('counter=%04d/%4d t/tmax=%.4f %s\n',counter,nt,t/tmax,fname); eval(['print -djpeg -r200 ' fname]); end