{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 }