% % 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=100; nx=1000; x=linspace(-L/2,L/2,nx)'; dx=x(2)-x(1); fprintf('L=%g nx=%d dx=%g\n',L,nx,dx); % Potential well : area under well is -alpha. % A square(ish) well with between -a/2 and a/2 and with area % (about) -alpha under it alpha = 11; a = 10; V0 = -alpha/a; V = V0*exp(-(x/(a/2)).^10); fprintf('Integral under potential = %g\n',sum(V*dx)); % 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*10/sig: sig=3; psi0 = exp(-(x+L/4).^2/sig^2).*exp(i*10/sig*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)); % 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 which is here % (#)*m*sig^2/h : a typical time scale for width sig is given by % m*sig^2/h from dimenstional analysis. The number 8 is dependent on what % we want to see and the initial momentum, etc. fs=20; counter=0; nt=90; tmax=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 ',t/(m*sig^2/(hbar*2*pi))) 'mw^2/ h']) 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