Assignment 9 solutions

1. > k:=200;

   

   > m:=1e-26;

   

   > L:=1.5e-10;

   

   > h := 6.6261e-34;

   

   > hbar := h/(2*Pi);

   

   > V := x -> k/2*(x-L/2)^2;

   

2. > phi1 := x -> sin(Pi*x/L) + a2*sin(2*Pi*x/L) + a3*sin(3*Pi*x/L);

   

   > Kip1 := -hbar^2/(2*m)*int(phi1(x)*diff(phi1(x),x$2),x=0..L);

   

   > Vip1 := int(phi1(x)^2*V(x),x=0..L);

   

   > ip1 := int(phi1(x)^2,x=0..L);

   

   > Evar1 := (Kip1+Vip1)/ip1;

   

   > s1:=solve({diff(Evar1,a2)=0,diff(Evar1,a3)=0},{a2,a3});

   

   > evalf(subs(s1[1],Evar1));

   

   > evalf(subs(s1[2],Evar1));

   

   > evalf(subs(s1[3],Evar1));

   

   > evalf(subs(s1[4],Evar1));

   

   > evalf(subs(s1[5],Evar1));

   

   The best solution is obtained with the third set of parameters. As indicated in the problem sheet, this set has a very small value of .

3. > phi2 := x -> sin(Pi*x/L) + a3*sin(3*Pi*x/L);

   

   > Kip2 := int(phi2(x)*(-hbar^2)/(2*m)*diff(phi2(x),x$2),x=0..L);

   

   > Vip2 := int(phi2(x)^2*V(x),x=0..L);

   

   > ip2 := int(phi2(x)^2,x=0..L);

   

   > Evar2 := (Kip2+Vip2)/ip2;

   

   > s2:=solve(diff(Evar2,a3)=0,a3);

   

   > subs(a3=s2[1],Evar2);

   

   > subs(a3=s2[2],Evar2);

   

   > subs(a3=s2[3],Evar2);

   

   The best solution for this trial wavefunction is obtained with the second value of . The energy is slightly better than what we obtained previously, i.e. should have been exactly zero in our previous attempt.

4. > phi3 := x -> sin(Pi*x/L) + a3*sin(3*Pi*x/L) + a5*sin(5*Pi*x/L);

   

   > Kip3 := int(phi3(x)*(-hbar^2)/(2*m)*diff(phi3(x),x$2),x=0..L);

   

   > Vip3 := int(phi3(x)^2*V(x),x=0..L);}{%

   

   > ip3 := int(phi3(x)^2,x=0..L);

   

   > Evar3 := (Kip3+Vip3)/ip3;

   

   > s3:=solve({diff(Evar3,a3)=0,diff(Evar3,a5)=0},{a3,a5});

   

   > subs(s3[1],Evar3);

   

   > subs(s3[2],Evar3);

   

   > subs(s3[3],Evar3);

   

   The lowest energy is obtained with the (1,3,5) wavefunction with the middle set of parameters.

5. > phi3n := x -> subs(s3[2],phi3(x)/sqrt(ip3));

   

   Verification:

   > int(phi3n(x)^2,x=0..L);

   

   > plot(phi3n,0..L);

   
6. Particle-in-a-box wavefunction:

   > plot(sqrt(2/L)*sin(Pi*x/L),x=0..L);

   

   Harmonic oscillator wavefunction:

   > alpha:=sqrt(k*m)/hbar;

   

   > plot((alpha/Pi)^(1/4)*exp(-alpha*(x-L/2)^2/2),x=0..L);

   

   Our wavefunction looks a little more like a harmonic oscillator wavefunction than like a particle-in-a-box wavefunction: There is relatively little density near the edges of the box and the main peak of the wavefunction is sharper. It would be interesting to pursue this calculation with more Fourier terms to see if the wavefunction more nearly approaches that of a harmonic oscillator. I expect it would: The fact that the harmonic oscillator wavefunction is so sharply peaked suggests to me that there is in fact very little interaction with the edges of the box and that the solution is in fact almost exactly the same as for a harmonic oscillator.

Bonus

Since the first excited state should have an antisymmetric wavefunction, we use a wavefunction which keeps only even-numbered terms of the Fourier series. Since the question isn't specific about how many terms to keep, let's just use two:

> phi4 := x -> sin(2*Pi*x/L) + a4*sin(4*Pi*x/L);

> Kip4 := int(phi4(x)*(-hbar^2)/(2*m)*diff(phi4(x),x$2),x=0..L);

> Vip4 := int(phi4(x)^2*V(x),x=0..L);

> ip4 := int(phi4(x)^2,x=0..L);

> Evar4 := (Kip4+Vip4)/ip4;

> s4:=solve(diff(Evar4,a4)=0,a4);

> subs(a4=s4[1],Evar4);

> subs(a4=s4[2],Evar4);

> phi4n := x -> subs(a4=s4[1],phi4(x)/sqrt(ip4));

Verification:

> int(phi4n(x)^2,x=0..L);

> plot(phi4n,0..L);