{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 283 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 286 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
287 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 288 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{CSTYLE "" -1 291 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 
292 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 293 "" 0 1 0 0 
0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 0 0 1 0 0 0 
0 0 0 0 }{CSTYLE "" -1 295 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 
2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 
{CSTYLE "" -1 -1 "Times" 1 16 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 
0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT -1 46 "The matrix representatio
n of quantum mechanics" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 
0 "" {TEXT -1 117 "We use the matrix representation of quantum mechani
cs to approximate the energy spectrum of an anharmonic oscillator." }}
{PARA 0 "" 0 "" {TEXT -1 119 "We make use of the particle-in-a-box bas
is. The box size has to be chosen sufficiently large so that the trunc
ation of " }{TEXT 264 1 "x" }{TEXT -1 39 "-space to a finite domain is
 justified." }}{PARA 0 "" 0 "" {TEXT -1 112 "An alternative would have
 been to choose a harmonic oscillator basis (as done in Quantum Mechan
ics using Maple)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 206 "We construct the eigenfunctions using the new LinearAlge
bra package of Maple 6. At the end we compare the approximate eigenfun
ction with a numerical solution to the differential-equation eigenvalu
e problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 279 "We solve the problem by calculating matrix elements in the coo
rdinate representation. Note, however, that a worksheet exists to carr
y out the calculations based on the commutation relations of quantum m
echanics. These worksheets are commut1.mws (Maple5), and define.mws (M
aple6)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
185 "We start with the basis. The box size will put limits on the accu
racy of the higher-lying eigenvalues/eigenfunctions. We also choose a \+
truncation size for the matrix eigenvalue problem." }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 41 "restart; Digits:=15: with(LinearAlgebra):" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "X:=8; N:=20;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "uB:=n->if is(n,even) then sqrt(1/X)
*sin(n*Pi*x/(2*X)) else sqrt(1/X)*cos(n*Pi*x/(2*X)) end if:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "plot([uB(1),uB(2),uB(3),uB(4
)],x=-X..X,color=[red,blue,green,black]);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 21 "int(uB(2)^2,x=-X..X);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 87 "The basis states are normalized; they represent the modes
 that the box can accomodate (" }{TEXT 258 1 "n" }{TEXT -1 86 "-1 coun
ts the number of nodes inside the box, and the waves vanish at the box
 edges (+" }{TEXT 257 1 "X" }{TEXT -1 6 " and -" }{TEXT 256 1 "X" }
{TEXT -1 280 "). What are the eigenenergies? We have kinetic energy on
ly, and evaluate it only over the allowed range; the result agrees wit
h the textbook formula. The kinetic energy matrix will be diagonal in \+
our basis, as the basis functions are eigenfunctions of the kinetic en
ergy operator." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "Tkin:=psi
->int(psi*(-1/2*diff(psi,x$2)),x=-X..X);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 50 "[Tkin(uB(1)),Tkin(uB(2)),Tkin(uB(3)),Tkin(uB(4))];" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "seq(n^2*Pi^2/8/X^2,n=1..4)
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "Now let us define our potent
ial of interest:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "V:=x^4;
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 170 "We define matrix elements wh
ere, again, space is restricted to a finite interval. This is imposed \+
by our choice of basis! We allow two different functions in bra and ke
t." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "Vpot:=(phi,psi)->int(
expand(phi*psi*V),x=-X..X);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "We
 are ready to assemble the hamiltonian matrix:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 14 "HM:=Matrix(N):" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 114 "We make use of the symmetry of the hamiltonian matrix: i
t allows to save almost a factor of 2 in computation time." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 126 "for i from 1 to N do: for j from 1
 to i do: HM[i,j]:=Vpot(uB(i),uB(j)); if i=j then HM[i,j]:=HM[i,j]+ Tk
in(uB(i)); fi: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "f
or i from 1 to N do: for j from i+1 to N do: HM[i,j]:=HM[j,i]: od: od:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "print(SubMatrix(HM,1..4
,1..4));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 251 "All integrals were c
alculated in closed form. The matrix has a sparse structure: the chose
n potential respects parity symmetry, and the basis consists of even/o
dd functions. Therefore, one could solve the problem of the even and o
dd eigenfunctions of " }{TEXT 259 1 "V" }{TEXT -1 100 " separately. We
 do keep it general, so that one may insert a non-symmetric potential \+
function above." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 126 "We cannot expect the matrix diagonalization to be carrie
d out symbolically. Therefore, we switch to floating-point evaluation:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "HMf:=map(evalf,HM):" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "evals:=Eigenvalues(HMf, output='lis
t'):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "ev_s:=sort(map(Re,evals));
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 260 11 "Exercise 1:" }}{PARA 0 "" 0 
"" {TEXT -1 173 "Carry out matrix diagonalizations with submatrices of
 the N-by-N Hamiltonian matrix, and observe the behaviour of the low-l
ying eigenvalues as a function of truncation size." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 123 "Now the eigenfun
ctions. We need a sorting procedure to arrange the result from the eig
envector calculation in proper order." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 36 "VE:=Eigenvectors(HMf,output='list'):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "Vp:=[seq([Re(VE[i][1]),VE[i][2],map
(Re,VE[i][3])],i=1..nops(VE))]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 344 "Min:=proc(x,y); if type(x,numeric) and type(y,numeric) then i
f x<=y then RETURN(true): else RETURN(false): fi; elif type(x,list) an
d type(y,list) and type(x[1],numeric) and type(y[1],numeric) then if x
[1]<=y[1] then RETURN(true): else RETURN(false): fi; elif convert(x,st
ring)<=convert(y,string) then RETURN(true): else RETURN(false): fi: en
d:\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "VEs:=sort(Vp,Min):
\n" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 252 "Suppose that we would like
 to see the eigenfunctions corresponding to the four lowest-lying eige
nvalues. The eigenvector for a given eigenvalue contains the expansion
 coefficients for the expansion of the eigenstate in terms of the chos
en basis states." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i f
rom 1 to 4 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "psi0:=add(VEs[i][
3][j]*uB(j),j=1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "No:=1/sqrt(
int(expand(psi0^2),x=-X..X));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "ph
i_a[i]:=add(No*VEs[i][3][j]*uB(j),j=1..N): od:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 66 "plot([seq(phi_a[i],i=1..4)],x=-5..5,color=[red
,blue,green,black]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "The four \+
functions should have 0,1,2,3 nodes. The additional nodes at large " }
{TEXT 261 1 "x" }{TEXT -1 44 " are clearly an artefact of the calculat
ion." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 73 "Let us check the accuracy of these approximate eigenfunctions o
ne by one:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "ev_s[1];" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "We can construct a shooting-method
 solution: we will use " }{TEXT 19 6 "dsolve" }{TEXT -1 60 " in numeri
c mode to integrate out the Schroedinger equation." }}{PARA 0 "" 0 "" 
{TEXT -1 59 "For the even-parity state we choose boundary conditions a
s " }{TEXT 19 17 "u(0)=1, D(u)(0)=0" }{TEXT -1 32 "; and for the odd-p
arity states " }{TEXT 19 17 "u(0)=0, D(u)(0)=1" }}{PARA 0 "" 0 "" 
{TEXT -1 128 "We use the energy as a trial value to ensure that the wa
vefunction vanishes at the 'boundary' (where the potential grows large
)." }}{PARA 0 "" 0 "" {TEXT -1 153 "This 'boundary' value is state-dep
endent. One iterates the procedure by hand. One can also write a bisec
tion algorithm to automate the eigenvalue search." }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 21 "IC:=u(0)=1,D(u)(0)=0;" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 11 "E_t:=0.668;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
39 "SE:=-1/2*diff(u(x),x$2)+(V-E_t)*u(x)=0;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 55 "sol:=dsolve(\{SE,IC\},u(x),numeric,output=listprocedu
re):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "u_s:=subs(sol,u(x)):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "P1:=plot('u_s(x)',x=0..2.5):
 plots[display](P1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "Now we ch
ange the normalization on the result of the matrix diagonalization to \+
make the comparison:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "P2:
=plot(phi_a[1]/subs(x=0,phi_a[1]),x=0..4,color=blue):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "plots[display](P1,P2);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT 262 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 
192 "Check the accuracy of the eigenenergies of the first and second e
xcited states, and compare the graphs of the eigenfunctions as obtaine
d by matrix diagonalization versus numerical integration." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 263 11 "Exercise 3:" }}
{PARA 0 "" 0 "" {TEXT -1 102 "Check the results obtained for different
 choices of the box size that defines the free-particle basis." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 478 "
It is possible to assess the accuracy of the matrix diagonalization wi
thout referring to the numerical calculation. One can calculate the ex
pectation value of the Hamiltonian in the coordinate representation fo
r the given eigenstate. This agrees with the eigenvalue obtained from \+
the matrix diagonalization. Then, one can ask about the distribution o
f energies in the state. The deviation from the average is calculated \+
from the matrix elements of the square of the Hamiltonian." }}{PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 209 "The functions to
 calculate these are defined below. However, they work only when the t
runcation size of the matrix is relatively small. Otherwise too much m
emory is consumed in the evaluation of the integrals." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "Havg:=psi->int(simplify(psi*(-1/2*d
iff(psi,x$2)+V*psi)),x=-X..X);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 23 "E1_avg:=Havg(phi_a[1]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 63 "H2me:=psi->int(simplify((-1/2*diff(psi,x$2)+V*psi)^2),x=-X..X)
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "Edev:=sqrt(abs(E1_avg^
2-H2me(phi_a[1])));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 283 "The devia
tion is much larger than what we would have expected from our comparis
on with the numerical value. Nevertheless, the number is proportional \+
to the accuracy for the given state: when the matrix size is increased
 the deviation from the average of the energy decreases for all " }
{TEXT 280 1 "x" }{TEXT -1 49 ". The large deviation is likely to come \+
from the " }{TEXT 281 1 "x" }{TEXT -1 124 "-range where the approximat
e wave function is oscillating. This can be verified by restricting th
e calculation to a smaller " }{TEXT 282 1 "x" }{TEXT -1 61 "-range whe
re the true ground-state eigenfunction is non-zero." }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 152 "Another way to explore
 the accuracy of the matrix diagonalization result is to graph the loc
al energy in position space. For an eigenstate the function " }{TEXT 
284 1 "f" }{TEXT -1 1 "(" }{TEXT 283 1 "x" }{TEXT -1 57 ") = (H*psi)/p
si should be equal to the eigenvalue at all " }{TEXT 285 1 "x" }{TEXT 
-1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "f:=psi->(-1/2*diff(psi,x$2))/psi+V;" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 42 "plot(f(phi_a[1]),x=0..3,view=[0..3,0..2])
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 451 "This graph should have a so
bering effect for any enthusiasm that may have developed after obtaini
ng a reasonably accurate eigenvalue, and an approximate eigenfunction \+
that appeared to follow the numerically calculated one: The local ener
gy calculation is a highly sensitive quantity. In particular, we notic
e that the basis-state representation in the 'particle-in-a-box' basis
 has difficulties with obtaining the correct fall-off of the wavefunct
ion." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 444 "
We should not be deterred too much by this mixture of success and disa
ppointment: we should be aware of the fact that eigenenergies are much
 easier to obtain than accurate eigenfunctions, and proceed with cauti
ously. Matrix diagonalization is often the only tool available to solv
e the Schroedinger equation. Usually we also have the possibility to i
mprove matters by adjusting the basis to the problem at hand, and ther
eby reducing the errors." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT 286 11 "Exercise 4:" }}{PARA 0 "" 0 "" {TEXT -1 218 "Grap
h the local energy for matrix diagonalization solutions differing by m
atrix size, and observe how the better converged calculations approxim
ate the eigenvalue locally. Repeat the calculation for a smaller value
 of " }{TEXT 287 1 "X" }{TEXT -1 7 " (with " }{TEXT 288 1 "X" }{TEXT 
-1 112 " larger than the point at which the numerically obtained wavef
unction goes to zero), and make your observations." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "We carry out the \+
calculation for the first excited state:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 8 "ev_s[2];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 
"IC:=u(0)=0,D(u)(0)=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "E
_t:=2.392;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "SE:=-1/2*diff(u(x),x$
2)+(V-E_t)*u(x)=0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "sol:=dsolve(
\{SE,IC\},u(x),numeric,output=listprocedure):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "u_s:=subs(sol,u(x)):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 68 "P1:=plot('u_s(x)',x=0..2.5): plots[display](P1,view=[
0..2.5,-1..1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 122 "We observe th
at the numerically obtained eigenvalues are slightly below the ones ob
tained from the matrix diagonalization." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 64 "P2:=plot(phi_a[2]/subs(x=0,diff(phi_a[2],x)),x=0..4,c
olor=blue):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "P3:=plot(V-E
_t,x=0..3,color=green):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "
plots[display](P1,P2,P3,view=[0..3,-2.5..1]);" }}}{EXCHG {PARA 0 "" 0 
"" {TEXT -1 124 "We have shifted the graph of the eigenfunction to ind
icate the location of the classical turning point (the crossing of the
 " }{TEXT 265 1 "x" }{TEXT -1 168 "-axis with the potential curve). Th
e plot shows size of the tunneling region, as well as the fact that th
e region where the numerical eigenfunction touches down to the " }
{TEXT 266 1 "x" }{TEXT -1 293 "-axis (it diverges for large x, due to \+
the finite accuracy of the trial eigenenergy, and due to numerical pro
blems) is well inside the classically forbidden region for the given s
tate. For higher-lying states one needs to go to larger x-values for t
he 'touchdown' point. Eventually, for large " }{TEXT 267 1 "n" }{TEXT 
-1 86 " the basis set for the matrix diagonalization will be inappropr
iate due to the finite " }{TEXT 268 1 "X" }{TEXT -1 7 "-value." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 714 "
It is important to explore the accuracy of states higher than the grou
nd states of the even- and odd-parity sectors. These states are harder
 to calculate, since one needs to find the locations of nodes that are
 not pre-determined. The eigenstates calculated from the matrix diagon
alization method are all mutually orthogonal. Given their approximate \+
nature, however, one cannot claim that the second excited state is ort
hogonal to the exact ground state. Experience shows that the accuracy \+
of the higher-lying eigenvalues (and their eigenfunctions) is substant
ially less than what was found for the two lowest states (for the same
 hamiltonian matrix). One can try to beat the problem by increasing th
e matrix size." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 92 "We demonstrate the problem on the example of the first ex
citation in the even-parity sector:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 8 "ev_s[3];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 
"IC:=u(0)=1,D(u)(0)=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "E
_t:=4.697;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "SE:=-1/2*diff(u(x),x$
2)+(V-E_t)*u(x)=0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "sol:=dsolve(
\{SE,IC\},u(x),numeric,output=listprocedure):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "u_s:=subs(sol,u(x)):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 66 "P1:=plot('u_s(x)',x=0..3): plots[display](P1,view=[0.
.3,-1.1..1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 122 "We observe that
 the numerically obtained eigenvalues are slightly below the ones obta
ined from the matrix diagonalization." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 56 "P2:=plot(phi_a[3]/subs(x=0,phi_a[3]),x=0..4,color=blu
e):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "P3:=plot(V-E_t,x=0..
3,color=green):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "plots[di
splay](P1,P2,P3,view=[0..3,-5..1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 549 "We observe that the third eigenvalue is not approximated well \+
by the 20-by-20 matrix diagonalization. This could be improved upon by
 simply choosing a smaller box. The numerical results for the wavefunc
tion combined with a demonstration of the vicinity of the classical tu
rning point (beyond which the tail of the wavefunction indicates that \+
tunneling takes place) indicate that the box size can be chosen to be \+
substantially smaller. This would lead to more accurate results, as th
e basis functions would be more flexible in the region of interest." }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT 295 96 "Solution of the numerical problem by bisection for the
 general case of non-symmetric potentials." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 "We start the integration from \+
" }{TEXT 272 1 "x" }{TEXT -1 3 "= -" }{TEXT 271 1 "s" }{TEXT -1 5 " an
d " }{TEXT 270 2 "x " }{TEXT -1 2 "= " }{TEXT 269 1 "s" }{TEXT -1 108 
", and try to match them at the origin, by requesting that u'(0)/u(0) \+
match from the left and from the right." }}{PARA 0 "" 0 "" {TEXT -1 
418 "The boundary condition at these two points distinguishes even- an
d odd-parity states. For even-parity states the derivatives of the wav
efunction are opposite to each other, for odd-parity states they are i
dentical. The choices determine the normalization of the propagated so
lution. Given that we match  u'(0)/u(0) this should not matter at all.
 Let us illustrate what happens as we integrate from the outside inwar
ds." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "s:=3.;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "eta:=3*10^(-4):" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 52 "IC1:=u(s)=0,D(u)(s)=eta: IC2:=u(-s)=0,D(u)(
-s)=-eta:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "E_t:=4.697;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 39 "SE:=-1/2*diff(u(x),x$2)+(V-E_t)*u(x)=0;" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 79 "sol1:=dsolve(\{SE,IC1\},u(x),numeric,output=l
istprocedure): u_1:=subs(sol1,u(x)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 79 "sol2:=dsolve(\{SE,IC2\},u(x),numeric,output=listprocedure): u_2:
=subs(sol2,u(x)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "P1:=pl
ot('u_1(x)',x=-3..3): plots[display](P1,view=[-3..3,-1..1]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "A matching method relies on the fa
ct that we wish the following properties for an eigenfunction:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 45 "normaliza
bility, i.e., vanishing function at " }{TEXT 274 1 "s" }{TEXT -1 6 " a
nd -" }{TEXT 273 1 "s" }{TEXT -1 1 ";" }}{PARA 0 "" 0 "" {TEXT -1 44 "
continuous and differentiable solution for -" }{TEXT 277 1 "s" }{TEXT 
-1 3 " < " }{TEXT 276 1 "x" }{TEXT -1 3 " < " }{TEXT 275 1 "s" }{TEXT 
-1 1 ";" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
199 "As the Schroedinger equation is of second order, we can produce t
wo linearly independent solutions that both depend on the same trial v
alue for the eigenenergy. One depends on initial conditions at -" }
{TEXT 279 1 "s" }{TEXT -1 40 ", the other on two intial conditions at \+
" }{TEXT 278 1 "s" }{TEXT -1 267 ". We calculate the difference betwee
n u'(0)/u(0) for the integrations from the left and from the right. We
 claim that this parameter crosses zero at an energy eigenvalue, and t
herefore we can invoke a root search. The latter is carried out by the
 bisection algorithm." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "Ma
tch:=proc(E_t) local SE,eta,IC1,IC2,sol1,sol2,u_1,up_1,u_2,up_2; globa
l V,s;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "SE:=-1/2*diff(u(x),x$2)+(
V-E_t)*u(x)=0; eta:=3*10^(-4):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "I
C1:=u(s)=0,D(u)(s)=eta: IC2:=u(-s)=0,D(u)(-s)=-eta:" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 79 "sol1:=dsolve(\{SE,IC1\},u(x),numeric,output=listpr
ocedure): u_1:=subs(sol1,u(x)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "
sol2:=dsolve(\{SE,IC2\},u(x),numeric,output=listprocedure): u_2:=subs(
sol2,u(x)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "up_1:=subs(sol1,diff
(u(x),x)): up_2:=subs(sol2,diff(u(x),x)):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "up_1(0)/u_1(0)-up_2(0)/u_2(0); end:" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 119 "The bisection algorithm is written with two pr
ocedures: a basic bisection step, and a driver which reports on progre
ss:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 115 "Bis1:=proc(x1,x2,x3
,f1,f2,f3);\nif evalf(f1*f3) < 0 then\nRETURN([x1,x3,f1,f3]);\nelse\nR
ETURN([x3,x2,f3,f2]); fi;\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 482 "Bisect:=proc(f,a,b) local res,x1,x2,x3,f1,f2,f3,i;x1:=a: x2:=
b: f1:=evalf(f(x1)): f2:=evalf(f(x2)):\nif evalf(f1*f2)>0 then RETURN(
\"No bracketed root\",f1,f2);\nelse\nx3:=0.5*(x1+x2); f3:=evalf(f(x3))
; fi;\nfor i from 1 to 50 do:\nres:=Bis1(x1,x2,x3,f1,f2,f3);\nx1:=res[
1]; f1:=res[3]; x2:=res[2]; f2:=res[4]; x3:=0.5*(x1+x2);\nif abs(x1-x2
) < 10^(-2) then print(\"Reached level 2\",x3); fi; if abs(x1-x2) < 10
^(-7) then RETURN(x3) fi;\nf3:=evalf(f(x3)); od:\nRETURN(\"Loop exhaus
ted\", x3); end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "s:=3;" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "Bisect(Match,4.65,4.75);" 
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 68 "Let us pick a non-symmetric pot
ential and calculate the eigenvalues:" }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "V:=x^4-2*x^3-x^2/2+2/2*
x+1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "P0:=plot(V,x=-5..5,
view=[-2..3,-1..5],color=green): plots[display](P0);" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 5 "s:=3;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "Bisect(Match,0,1);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 6 "E1:=%;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 109 "To gra
ph the solution we use the obtained best eigenvalue, and integrate the
 Schroedinger equation once more:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 15 "eta:=1*10^(-6):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 51 "IC1:=u(s)=0,D(u)(s)=eta: IC2:=u(-s)=0,D(u)(-s)=eta:" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "E_t:=E1;" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 39 "SE:=-1/2*diff(u(x),x$2)+(V-E_t)*u(x)=0;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "sol1:=dsolve(\{SE,IC1\},u(x),numeri
c,output=listprocedure): u_1:=subs(sol1,u(x)):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 79 "sol2:=dsolve(\{SE,IC2\},u(x),numeric,output=listproce
dure): u_2:=subs(sol2,u(x)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 77 "P1:=plot('u_1(x)',x=-s..s,color=red): P2:=plot('u_2(x)',x=-s..s,
color=blue): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "P00:=plot(
E1,x=-3..3,color=black):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 48 
"plots[display](P00,P0,P1,P2,view=[-3..3,-1..6]);" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 233 "We have chosen the initial condition for solution
 2 to be such that a positive ground-state wavefunction is calculated.
 Solution 1 (which is integrated left-to-right), is small, and cannot \+
be seen on this scale (it is covered by the " }{TEXT 289 1 "x" }{TEXT 
-1 60 "-axis). Both solutions are normalized through the choice of " }
{TEXT 19 3 "eta" }{TEXT -1 135 ". It is interesting to observe how the
 wavefunction reaches through the classically forbidden regime to exte
nd over the second minimum." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT 290 11 "Exercise 5:" }}{PARA 0 "" 0 "" {TEXT -1 115 "C
alculate some of the excited eigenstates for this potential by supplyi
ng other bracketing values to the procedure " }{TEXT 19 6 "Bisect" }
{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT 291 11 "Exercise 6:" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
167 "Modify the potential function such that there be one broad, and o
ne sharp minimum with approximately equal depth by choosing different \+
coeffients for the monomials in " }{TEXT 292 1 "x" }{TEXT -1 162 ". Ca
lculate the ground state and observe its shape. When does it become do
uble-humped, and when does the second hump (over the narrower potentia
l well) disappear?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT 293 23 "Matrix diagonalization:" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 172 "We can get an idea about
 the location of the low-lying eigenvalues by carrying out the matrix \+
diagonalization. We simply copy the lines from the beginning of the wo
rksheet." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "HM:=Matrix(N):
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 114 "We make use of the symmetry \+
of the hamiltonian matrix: it allows to save almost a factor of 2 in c
omputation time." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 126 "for i \+
from 1 to N do: for j from 1 to i do: HM[i,j]:=Vpot(uB(i),uB(j)); if i
=j then HM[i,j]:=HM[i,j]+ Tkin(uB(i)); fi: od: od:" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 71 "for i from 1 to N do: for j from i+1 to N d
o: HM[i,j]:=HM[j,i]: od: od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 31 "print(SubMatrix(HM,1..4,1..4));" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 78 "The Hamiltonian matrix is no longer sparse, as the potent
ial is not symmetric." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "HM
f:=map(evalf,HM):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "evals:=Eigenva
lues(HMf, output='list'):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "ev_s:=
sort(map(Re,evals));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "VE:
=Eigenvectors(HMf,output='list'):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 66 "Vp:=[seq([Re(VE[i][1]),VE[i][2],map(Re,VE[i][3])],i=1
..nops(VE))]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "VEs:=sort(
Vp,Min):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to
 4 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "psi0:=add(VEs[i][3][j]*uB
(j),j=1..N):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "No:=1/sqrt(int(expa
nd(psi0^2),x=-X..X));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "phi_a[i]:=
add(No*VEs[i][3][j]*uB(j),j=1..N): od:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 66 "plot([seq(phi_a[i],i=1..4)],x=-5..5,color=[red,blue,g
reen,black]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 138 "We observe that
 the matrix diagonalization gets the right idea about the ground-state
 eigenfunction (the eigenvalue is off by nearly 3 %)." }}{PARA 0 "" 0 
"" {TEXT -1 83 "The first excited state has a bigger hump over the sha
llower part of the potential." }}{PARA 0 "" 0 "" {TEXT -1 180 "The hig
her states are not unlike the quartic harmonic oscillator eigenstates,
 but shifted to the right. For higher eigenenergies the particles expl
ore mostly the steep rise in the " }{TEXT 19 3 "x^4" }{TEXT -1 23 "-pa
rt of the potential." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT 294 11 "Exercise 7:" }}{PARA 0 "" 0 "" {TEXT -1 181 "Choose a \+
potential with a deep and narrow potential well, and make a detailed c
omparison between the ground state as obtained by matrix-diagonalizati
on, and by the numerical method." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 0 "" }}}}{MARK "115" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 
}{PAGENUMBERS 0 1 2 33 1 1 }