- We will need the particle-in-a-box wavefunctions:
> psi1d := (n,x) -> sqrt(2/L)*sin(n*Pi*x/L);
The probability that the particle will be found between L/3
and L/2 for n=1 is
> evalf(int(psi1d(1,x)^2,x=L/3..L/2));
For n=2, we get
> evalf(int(psi1d(2,x)^2,x=L/3..L/2));
For n=1, the probability density has a maximum at
x=L/2 and rises monotonically in the interval
studied. On the other hand, the probability density decreases
monotonically to a node at x=L/2 in this interval
for n=3. -
The wavefunctions are orthogonal if their inner product is zero:
> int(psi1d(3,x)*psi1d(212,x),x=0..L);
Wavefunctions corresponding to different energies are always orthogonal,
so we really didn't need to do this calculation. -
The potential energy is treated as a perturbation to the particle-in-a-box
solutions.
> V := x -> h^2/(8*m*L^4)*(x-L/2)^2;
The energy correction for n=1 is
> deltaE1 := int(psi1d(1,x)*V(x)*psi1d(1,x),x=0..L);
so the total energy of the perturbed system for n=1 is
> E1 := h^2/(8*m*L^2) + deltaE1;
For n=2, we get
> deltaE2 := int(psi1d(2,x)*V(x)*psi1d(2,x),x=0..L);
> E2 := 2^2*h^2/(8*m*L^2) + deltaE2;
In the unperturbed system, the energy difference between these two
states is 
:
> evalf(3*h^2/(8*m*L^2));
In the perturbed system, we have
> evalf(E2-E1);
This is bigger than the energy gap in the unperturbed system. -
> psi2d := (nx,ny,x,y) -> 2/L*sin(nx*Pi*x/L)*sin(ny*Pi*y/L);
-
We need to integrate the wavefunction over the square. For
nx=1, ny=2, we get
> evalf(int(int(psi2d(1,2,x,y)^2,x=0..L/3),y=L/3..2*L/3));
and for nx=2, ny=1
> evalf(int(int(psi2d(2,1,x,y)^2,x=0..L/3),y=L/3..2*L/3));
-
> int(int(psi2d(1,2,x,y)*x^2*y^2*psi2d(1,2,x,y),x=0..L),y=0..L);
