a r X i v :0704.1622v 1 [p h y s i c s .e d -p h ] 12 A p r 2007
MATLAB codes for teaching quantum physics:Part 1
R.Garcia,∗ A.Zozulya,and J.Stickney
Department of Physics,Worcester Polytechnic Institute,Worcester,MA 01609
(Dated:February 1,2008)
Among the ideas to be conveyed to students in an introductory quantum cour,we have the pivotal idea championed by Dirac that functions correspond to column vectors (kets)and that differential operators correspond to matrices (ket-bras)acting on tho vectors.The MATLAB (matrix-laboratory)programming environment is especially uful in conveying the concepts to students becau it is geared towards the type of matrix manipulations uful in solving introductory quantum physics problems.In this article,we share MATLAB codes which have been developed at WPI,focusing on 1D problems,to be ud in conjunction with Griffiths’introductory text.
Two key concepts underpinning quantum physics are the Schrodinger equation and the Born probability equa-tion.In 1930Dirac introduced bra-ket notation for state vectors and operators.1This not
ation emphasized and clarified the role of inner products and linear function spaces in the two equations and is fundamental to our modern understanding of quantum mechanics.The Schrodinger equation tells us how the state Ψof a particle evolves in time.In bra-ket notation,it reads
i d
2m
∂ x
2
Ψ( x ,t )+U ( x )
x |(3)
acting on the space of kets.
While an expert will necessarily regard Eqs.(1-3)as a great simplification when thinking of the content of quan-tum physics,the novice often understandably reels under the weight of the immen
abstraction.We learn much about student thinking from from the answers given by our best students.For example,we find a common error
when studying 1D quantum mechanics is a student treat-ing Ψ(x )and |Ψ interchangeably,ignoring the fact that the first is a scalar but the ket corresponds to a column vector.For example,they may write incorrectly
p ||Ψ |x =|Ψ p |x (incorrect!)(4)
or some similar abberation.To avoid the types of mis-conceptions,a number of educators and textbook authors have stresd incorporating a numerical calculation as-pect to quantum cours.4,5,6,7,8,9The motive is simple.Anyone who has done numerical calculations can’t help but regard a ket |Ψ as a column vector,a bra Ψ|as a row vector and an operator H as a matrix becau that is how they concretely reprented in the computer.Intro-ducing a computational aspect to the cour provides one further benefit:it gives the beginning quantum student the n that he or she is being empowered to solve real problems that may not have simple,analytic solutions.With the motivations in mind,we have developed MATLAB codes 10for solving typical 1D problems found in the first part of a junior level quantum cour bad on Griffith’s book.11We cho MATLAB for our p
ro-gramming environment becau the MATLAB syntax is especially simple for the typical matrix operations ud in 1D quantum mechanics problems and becau of the ea of plotting functions.While some MATLAB numeri-cal recipes have previously been published by others,12,13the exercis we share here are special becau they em-phasize simplicity and quantum pedagogy,not numerical efficiency.Our point has been to provide exercis which show students how to numerically solve 1D problems in such a way that emphasizes the column vector aspect of kets,the row vector aspect of bras and the matrix aspect of operators.Exercis using more efficient MATLAB ODE solvers or finite-element techniques are omitted be-cau they do not rve this immediate purpo.In part II of this article,we hope to share MATLAB codes which can be ud in conjunction with teaching topics pertain-ing to angular momentum and non-commuting obrv-ables.
2
I.FUNCTIONS AS VECTORS
To start students thinking of functions as column vector-like objects,it is very uful to introduce them to plotting and integrating functions in the MATLAB environment.Interestingly enough,the plot comm
and in MATLAB takes vectors as its basic input element.As shown in Program1below,to plot a function f(x)in MATLAB,wefirst generate two vectors:a vector of x values and a vector of y values where y=f(x).The com-mand plot(x,y,′r′)then generates a plot window con-taining the points(x i,y i)displayed as red points(′r′). Having specified both x and y,to evaluate the definite integral L−L ydx,we need only sum all the y values and multiply by dx.
%**************************************************************** %Program1:Numerical Integration and Plotting using MATLAB
%**************************************************************** N=1000000;%No.of points
猪湾L=500;%Range of x:from-L to L
x=linspace(-L,L,N)’;%Generate column vector with N
%x values ranging from-L to L
dx=x(2)-x(1);%Distance between adjacent points
%Alternative Trial functions:
%To lect one,take out the comment command%at the beginning. %y=exp(-x.^2/16);%Gaussian centered at x=0
%y=((2/pi)^0.5)*exp(-2*(x-1).^2);%Normed Gaussian at x=1
%y=(1-x.^2).^-1;%Symmetric fcn which blows up at x=1 %y=(cos(pi*x)).^2;%Cosine fcn
%y=exp(i*pi*x);%Complex exponential
%y=sin(pi*x/L).*cos(pi*x/L);%Product of sinx times cosx
%y=sin(pi*x/L).*sin(pi*x/L);%Product of sin(nx)times sin(mx)
%A=100;y=sin(x*A)./(pi*x);%Rep.of delta fcn
A=20;y=(sin(x*A).^2)./(pi*(A*x.^2));%Another rep.of delta fcn %Obrve:numerically a function y(x)is reprented by a vector!
%Plot a vector/function
plot(x,y);%Plots vector y vs.x
%plot(x,real(y),’r’,x,imag(y),’b’);%Plots real&imag y vs.x axis([-2207]);%Optimized axis parameters for sinx^2/pix^2 %axis([-22-840]);%Optimized axis parameters for sinx/pix
%Numerical Integration
sum(y)*dx%Simple numerical integral of y
trapz(y)*dx%Integration using trapezoidal rule
II.DIFFERENTIAL OPERATORS AS
MATRICES
Just as f(x)is reprented by a column vector|f in the computer,for numerical purpos a differential opera-torˆD acting on f(x)is rerented by a matrix D that acts on|f .As illustrated in Program2,MATLAB provides many uful,intuitive,well-documented commands for generating matrices D that correspond to a givenˆD.10 Two examples are the commands ones and diag.The command ones(a,b)generates an a×b matrix of ones. The command diag(A,n)generates a matrix with the el-ements of the vector A placed along the n th diagonal and zeros everywhere el.The central diagonal corresponds to n=0,the diagonal above the center one corresponds to n=).
An exerci we suggest is for students to verify that the derivative matrix is not Hermitian while the deriva-tive matrix times the imaginary number i is.This can be very valuable for promoting student understanding if done in conjunction with the proof usually given for the differential operator.
%**************************************************************** %Program2:Calculate first and cond derivative numerically
%showing how to write differential operator as a matrix
%**************************************************************** %Parameters for solving problem in the interval0<x<L
L=2*pi;%Interval Length
N=100;%No.of coordinate points
x=linspace(0,L,N)’;%Coordinate vector
dx=x(2)-x(1);%Coordinate step
酒店前台实习周记
%Two-point finite-difference reprentation of Derivative
D=(diag(ones((N-1),1),1)-diag(ones((N-1),1),-1))/(2*dx);
%Next modify D so that it is consistent with f(0)=f(L)=0
D(1,1)=0;D(1,2)=0;D(2,1)=0;%So that f(0)=0
D(N,N-1)=0;D(N-1,N)=0;D(N,N)=0;%So that f(L)=0
%Three-point finite-difference reprentation of Laplacian
Lap=(-2*diag(ones(N,1),0)+diag(ones((N-1),1),1)...
+diag(ones((N-1),1),-1))/(dx^2);
%Next modify Lap so that it is consistent with f(0)=f(L)=0 Lap(1,1)=0;Lap(1,2)=0;Lap(2,1)=0;%So that f(0)=0 Lap(N,N-1)=0;Lap(N-1,N)=0;Lap(N,N)=0;%So that f(L)=0
%To verify that D*f corresponds to taking the derivative of f
%and Lap*f corresponds to taking a cond derviative of f,
%define f=sin(x)or choo your own f
f=sin(x);
%And try the following:
Df=D*f;Lapf=Lap*f;
plot(x,f,’b’,x,Df,’r’,x,Lapf,’g’);
axis([05-1.11.1]);%Optimized axis parameters
%To display the matrix D on screen,simply type D D%Displays the matrix D in the workspace
Lap%Displays the matrix Lap
III.INFINITE SQUARE WELL
When solving Eq.(1),the method of paration of vari-ables entails that as an intermediate step we look for the parable solutions
|ΨE(t) =|ΨE(0) exp(−iEt/ )(5) where|ΨE(0) satisfies the time-independent Schrodinger equation
H|ΨE(0) =E|ΨE(0) .(6) In solving Eq.(6)we are solving for the eigenvalues E and eigenvectors|ΨE(0) of H.In MATLAB,the com-mand[V,E]=eig(H)does precily this:it generates two matrices.Thefirst matrix V has as its columns the eigenvectors|ΨE(0) .The cond matrix E is a diago-nal matrix with the eigenvalues E i corresponding to the eigenvectors|ΨE i(0) placed along the central diagonal. We can u the command E=diag(E)to convert this matrix into a column vector.In Program3,we solve for
3
the eigenfunctions and eigenvalues for the infinite square well Hamiltonian.For brevity,we omit the commands tting the parameters L,N,x,and dx.
%**************************************************************** %Program3:Matrix reprentation of differential operators,
%Solving for Eigenvectors&Eigenvalues of Infinite Square Well %**************************************************************** %For brevity we omit the commands tting the parameters L,N,
%x and dx;We also omit the commands defining the matrix Lap.
%The would be the same as in Program2above.
%Total Hamiltonian where hbar=1and m=1
hbar=1;m=1;H=-(1/2)*(hbar^2/m)*Lap;
%Solve for eigenvector matrix V and eigenvalue matrix E of H [V,E]=eig(H);
%Plot lowest3eigenfunctions
plot(x,V(:,3),’r’,x,V(:,4),’b’,x,V(:,5),’k’);shg;
E%display eigenvalue matrix
diag(E)%display a vector containing the eigenvalues
Note that in the MATLAB syntax the object V(:,3) specifies the column vector consisting of all the elements in column3of matrix V.Similarly V(2,:)is the row vector consisting of all elements in row2of V;V(3,1)is the element at row3,column1of V;and V(2,1:3)is a row vector consisting of elements V(2,1),V(2,2)and V(2,3).
IV.ARBITRARY POTENTIALS
Numerical solution of Eq.(1)is not limited to any par-ticular potential.Program4gives example MATLAB codes solving the time independent Schrodinger equa-tion forfinite square well potentials,the harmonic os-cillator potential and even for potentials that can only solved numerically such as the quartic potential U=x4. In order to minimize the amount of RAM required,the codes shown make u of spar matrices,where only the non-zero elements of the matrices are stored.The commands for spar matrices are very similar to tho for non-spar matrices.For example,the command [V,E]=eigs(H,nmodes..)provides the nmodes lowest energy eigenvectors V of of the spar matrix H.
Fig.1shows the plot obtained from Program4for the potential U=1
2 (7) with n=integer andω= m=10rad/s.
散步预习
V.A NOTE ON UNITS IN OUR PROGRAMS
When doing numerical calculations,it is important to minimize the effect of rounding errors by choosing units such that the variables ud in the simulation are of the %*******************************
********************************* %Program4:Find veral lowest eigenmodes V(x)and
%eigenenergies E of1D Schrodinger equation
%-1/2*hbar^2/m(d2/dx2)V(x)+U(x)V(x)=EV(x)
%for arbitrary potentials U(x)
%**************************************************************** %Parameters for solving problem in the interval-L<x<L
%PARAMETERS:
L=5;%Interval Length
N=1000;%No of points
x=linspace(-L,L,N)’;%Coordinate vector
dx=x(2)-x(1);%Coordinate step
%POTENTIAL,choo one or make your own
U=1/2*100*x.^(2);%quadratic harmonic oscillator potential %U=1/2*x.^(4);%quartic potential
%Finite square well of width2w and depth given
%w=L/50;
%U=-500*(heaviside(x+w)-heaviside(x-w));
%Two finite square wells of width2w and distance2a apart
%w=L/50;a=3*w;
%U=-200*(heaviside(x+w-a)-heaviside(x-w-a)...
%+heaviside(x+w+a)-heaviside(x-w+a));
%Three-point finite-difference reprentation of Laplacian
%using spar matrices,where you save memory by only
%storing non-zero matrix elements
e=ones(N,1);Lap=spdiags([e-2*e e],[-101],N,N)/dx^2;
%Total Hamiltonian
hbar=1;m=1;%constants for Hamiltonian
H=-1/2*(hbar^2/m)*Lap+spdiags(U,0,N,N);
%Find lowest nmodes eigenvectors and eigenvalues of spar matrix nmodes=3;options.disp=0;
[V,E]=eigs(H,nmodes,’sa’,options);%find eigs
[E,ind]=sort(diag(E));%convert E to vector and sort low to high V=V(:,ind);%rearrange corresponding eigenvectors
%Generate plot of lowest energy eigenvectors V(x)and U(x)
Usc=U*max(abs(V(:)))/max(abs(U));%rescale U for plotting plot(x,V,x,Usc,’--k’);%plot V(x)and rescaled U(x)
%Add legend showing Energy of plotted V(x)
lgnd_str=[repmat(’E=’,nmodes,1),num2str(E)];
legend(lgnd_str)%place lengend string on plot
shg
order of unity.In the programs prented here,our fo-cus being undergraduate physics students,we wanted to avoid unnecessarily complicating matters.To make the equations more familiar to the students,we explicitly left constants such as in the formulas and cho units such that =1and m=1.We recognize that others may have other opinions on how to address this issue.An al-ternative approach ud in rearch is to recast the equa-tions in terms of dimensionless variables,for example by rescaling the energy to make it dimensionless by express-ing it in terms of some characteristic energy in the prob-lem.In a more advanced cour for graduate students or in a cour in numerical methods,such is an approach which would be preferable.
VI.TIME DEPENDENT PHENOMENA The parable solutions|ΨE(t) are only a subt of all possible solutions of Eq.(1).Fortunately,they are complete t so that we can construct the general solution
4
FIG.1:Output of Program4,which plots the energy eigen-
functions V(x)a nd a scaled version of the potential U(x)=
1/2·100·x2.The corresponding energies displayed within
thefigure legend,4.99969,14.9984and24.9959,are,within
rounding error,precily tho expected from Eq.(7)for the
three lowest-energy modes.
小寒吃什么via the linear superposition
|Ψ(t) = E a E|ΨE(0) exp(−iEt/ )(8)
where a E are constants and the sum is over all possible
values of E.The important difference between the p-
arable solutions(5)and the general solution(8)is that
the probability densities derived from the general solu-
tions are time-dependent whereas tho derived from the
parable solutions are not.A very apt demonstration
of this is provided in the Program5which calculates the
time-dependent probability densityρ(x,t)for a particle
trapped in a a pair offinite-square wells who initial
state|Ψ(0) is t equal to the the equally-weighted su-
perposition
|Ψ(0) =12(|ΨE0 +|ΨE1 )(9)
of the ground state|ΨE
0 andfirst excited state|ΨE1
of the double well system.As snapshots of the program
output show in Fig.2,the particle is initially completely
localized in the rightmost well.However,due to E0=E1,
the probability density
1
ρ(x,t)=
5
FIG.2:The probability densityρ(x,t)output from Program 5for a particle trapped in a pair offinite-height square well potentials that are cloly adjacent.The initial state of the particle is chon to be|Ψ(0) ∼|E0 +|E1 .Shown isρ(x,t) plotted for times T={0,0.25,0.5,0.75,1.0}×2π /(E1−E2). VII.W A VEPACKETS AND STEP POTENTIALS Wavepackets are another time-dependent phenomenon encountered in undergraduate quantum mechanics for which numerical solution techniques have been typically advocated in the hopes of promoting intuitive acceptance and understanding of approximations necessarily invoked in more formal,analytic treatments.Program6calcu-lates and displays the time evolution of a wavepacket for one of two possible potentials,either U=0or a step potential U=U oΘ(x−L).The initial wavepacket is generated as the Fast Fourier Transform of a Gauss
ian momentum distribution centered on a particular value of the wavevector k o.Becau the wavepacket is com-pod of a distribution of different k s,different parts of the wavepacket move with different speeds,which leads to the wave packet spreading out in space as it moves. While there is a distribution of velocities within the wavepacket,two velocities in particular are uful in char-acterizing it.The pha velocity v w=ω/k=E/p= k o/2m is the velocity of the plane wave component which has wavevector k o.The group velocity v g= k o/m is the velocity with which the expectation value<x> moves and is the same as the classical particle velocity as-sociated with the momentum p= k.Choosing U=0, students can modify this program to plot<x>vs t. They can extract the group velocity from their numerical simulation and obrve that indeed v g=2v w for a typ-ical wave packet.Students can also obrve that,while v g matches the particle speed from classical mechanics, the wavepacket spreads out as time elaps.
In Program6,we propagate the wave function forward via the formal solution湿气过重
|Ψ(t) =exp(−iHt/ )|Ψ(0) ,(13) where the Hamiltonian matrix H is in the exponential. This solution is equivalent to Eq.(6),as as can be shown by simple substitution.Moreover,MATLAB has no trou-ble exponentiating the matrix that numerically repre-nting the Hamiltonian operator as long as the matrix is small enough tofit in the available computer memory. Even more interestingly,students can
u this method to investigate scattering of wavepackets from various po-tentials,including the step potential U=U oΘ(x−L/2). In Fig.3,we show the results of what happens as the wavepacket impinges on the potential barrier.The pa-rameters characterizing the initial wavepacket have been deliberately chon so that the wings do not fall outside the simulation area and initially also do not overlap the barrier on the right.If E ≪U o,the wavepacket is completely reflected from the barrier.If E ≈U o,a portion of the wave is is reflected and a portion is trans-mitted through.If E ≫U o,almost all of the wave is transmitted.
横看
In Fig.4we compare the reflection probability R cal-culated numerically using Program6with R calculated by averaging the single-mode11expression
R(E)= √E−U o
E+
√
6
%**************************************************************** %Program6:Wavepacket propagation usi
ng exponential of H
%**************************************************************** %Parameters for solving the problem in the interval0<x<L
L=100;%Interval Length
N=400;%No of points
x=linspace(0,L,N)’;%Coordinate vector
dx=x(2)-x(1);%Coordinate step
%Parameters for making intial momentum space wavefunction phi(k) ko=2;%Peak momentum
a=20;%Momentum width parameter
dk=2*pi/L;%Momentum step
km=N*dk;%Momentum limit
k=linspace(0,+km,N)’;%Momentum vector
%Make psi(x,0)from Gaussian kspace wavefunction phi(k)using
%fast fourier transform:
phi=exp(-a*(k-ko).^2).*exp(-i*6*k.^2);%unnormalized phi(k)
psi=ifft(phi);%multiplies phi by expikx and integrates vs.x psi=psi/sqrt(psi’*psi*dx);%normalize the psi(x,0) %Expectation value of for the parameters
%chon above<E>=2.062.
avgE=phi’*0.5*diag(k.^2,0)*phi*dk/(phi’*phi*dk);
%CHOOSE POTENTIAL U(X):Either U=0OR
%U=step potential of height avgE that is located at x=L/2
%U=0*heaviside(x-(L/2));%free particle wave packet evolution U=avgE*heaviside(x-(L/2));%scattering off step potential
%Finite-difference reprentation of Laplacian and Hamiltonian
e=ones(N,1);Lap=spdiags([e-2*e e],[-101],N,N)/dx^2;
H=-(1/2)*Lap+spdiags(U,0,N,N);
%Parameters for computing psi(x,T)at different times0<T<TF NT=200;%No.of time steps莲菜怎么凉拌
TF=29;T=linspace(0,TF,NT);%Time vector
dT=T(2)-T(1);%Time step
hbar=1;
%Time displacement operator E=exp(-iHdT/hbar)
E=expm(-i*full(H)*dT/hbar);%time desplacement operator
%*************************************************************** %Simulate rho(x,T)and plot for each T
%*************************************************************** for t=1:NT;%time index for loop
%calculate probability density rho(x,T)
psi=E*psi;%calculate new psi from old psi
rho=conj(psi).*psi;%rho(x,T)
plot(x,rho,’k’);%plot rho(x,T)vs.x
axis([0L00.15]);%t x,y axis parameters for plotting
xlabel(’x[m]’,’FontSize’,24);
ylabel(’probability density[1/m]’,’FontSize’,24);
pau(0.05);%pau between each frame displayed
end
%Calculate Reflection probability
R=0;
for a=1:N/2;
R=R+rho(a);
end
R=R*dx
numerical simulation for E /U o≈1.This discrepancy can be reduced significantly by increasing the number of points in the simulation to1600but only at the cost of significantly slowing down the speed of the computation. For our purpos,the importance comparing the analyt-ical and numerical calculations is that it gives student a baline from which to
长春几线城市
form an opinion or intuition regarding the accuracy of Eq.(14).FIG.3:Output of Program6showing a wavepacket encoun-tering step potential of height∼ E located at x/L=0.5at different times.
VIII.CONCLUSIONS
One benefit of incorporating numerical simulation into the teaching of quantum mechanics,as we have men-tioned,is the development of student intuition.Another is showing students that non-ideal,real-world problems can be solved using the concepts they learn in the class-room.However,our experimentation incorporating the simulations in quantum physics at WPI during the past year has shown us that the most important benefit is a type of side-effect to doing numerical simulation:the ac-ceptance on an intuitive level by the student that func-tions are like vectors and differential operators are like matrices.While in the prent paper,we have only had sufficient space to prent the most illustrative MATLAB
