Notebook to demonstrate the variational method for the infinite square well

This deals with the well in which the potential is zero between 0 and a, and infinite elsewhere.  The wave function must therefore vanish except in the region 0-a. Since the potential vanishes in this region, only the term (p^2)/2m  contributes to the energy.  In what follows we will use<Ψ| (p^2)|Ψ>= (ℏ^2)<((d Ψ)/(d x))|((d Ψ)/(d x))>.  The true energy values are (Subscript[E, n])=(n^2)(π^2) (ℏ^2)/2(ma^2) (with the ground state energy being (Subscript[E, 1])) and from now on we neglect the overall factor of  (ℏ^2)/2m.  Thus we should get (Subscript[E, 1])≥(π^2)/(a^2) and (Subscript[E, 2])≥4(π^2)/(a^2).

First approximation to the ground state

(Subscript[Ψ, 1])[x_]=x(a-x);

This trial wavefunction obeys the boundary conditions (Ψ[0]=Ψ[a]=0), has no nodes, and is symmetric about x=a - all attributes we expect of the true ground state.

Plot[(Subscript[Ψ, 1])[x]/.a->1,{x,0,1},ImageSize->200]

Simplify[Integrate[D[(Subscript[Ψ, 1])[x],x]^2,{x,0,a}]/Integrate[(Subscript[Ψ, 1])[x]^2,{x,0,a}]]

%/((π^2)/(a^2))//N

So this very simple trial wavefunction has given us a lower bound on the energy just over 1% above the true value.

Approximation to the first excited state

(Subscript[Ψ, 2])[x_]=x(a-x)(2x-a);

Here we have chosen a trial wave function with is antisymmetric about x=a/2

Plot[(Subscript[Ψ, 2])[x]/.a->1,{x,0,1},ImageSize->200]

Simplify[Integrate[D[(Subscript[Ψ, 2])[x],x]^2,{x,0,a}]/Integrate[(Subscript[Ψ, 2])[x]^2,{x,0,a}]]

%/(4(π^2)/(a^2))//N

Improved approximation to the ground state

(Subscript[Ψ, 3])[x_]=((x/a)(1-x/a)+b(x/a)^2(1-x/a)^2);

This is just the first function I thought of which continues to satisfy the constraints discussed above, but also has a free parameter b.  (Unless if b is too negative, in which case  the function will develop nodes)

Eb[b_]=Simplify[Integrate[D[(Subscript[Ψ, 3])[x],x]^2,{x,0,a}]/Integrate[(Subscript[Ψ, 3])[x]^2,{x,0,a}]]

Plot[a^2Eb[b],{b,-.2,3.5},ImageSize->250,AxesOrigin->{0,Pi^2}]

This isn' t an improvement for all values of b, but it is for those between (roughly) 0 and 3.

min=FullSimplify[Minimize[(a^2)Eb[b],b]]

min//N

min[[1]]/((π^2))//N

Good to one part in (10^5)!

nint=Integrate[(((Subscript[Ψ, 3])[x])^2),{x,0,a}]

Plot[{Sqrt[2] Sin[x Pi], Sqrt[30] x(1-x),(Subscript[Ψ, 3])[x]/Sqrt[nint]/.a->1/.min[[2]]},{x,0,1},AxesLabel->{x/a,Ψ},PlotStyle->{Black,Red,Blue},PlotLabel->"True and aproximate wave functions"]

Plot[{ Sqrt[30] x(1-x)-Sqrt[2] Sin[x Pi],((Subscript[Ψ, 3])[x]/Sqrt[nint]-Sqrt[2] Sin[x Pi])/.a->1/.min[[2]]},{x,0,1},AxesLabel->{x/a,"Ψ-(Subscript[Ψ, 0])"},PlotStyle->{Red,Blue},PlotLabel->"Deviations of approximate wavefunctions from the true one"]

We can hardly see the blue curve, so lets multiply it by 5:

Plot[{ Sqrt[30] x(1-x)-Sqrt[2] Sin[x Pi],5((Subscript[Ψ, 3])[x]/Sqrt[nint]-Sqrt[2] Sin[x Pi])/.a->1/.min[[2]]},{x,0,1},AxesLabel->{x/a,"Ψ-(Subscript[Ψ, 0])"},PlotStyle->{Red,Blue},PlotLabel->"Deviations of approximate wavefunctions from the true one"]

Challenge

Now you see if you can do better for (Subscript[E, 2])...

Higher excited states

What about (Subscript[E, 3])?   Let's find a wavefuction with the form of (Subscript[Ψ, 3])[x] but with a new value of b so it is orthogonal to the best trial wave function for the ground state:

Integrate[((Subscript[Ψ, 3])[x]/.{a->1,b->b2})((Subscript[Ψ, 3])[x]/.{a->1}),{x,0,1}]

Solve[%==0,b2]

FullSimplify[Flatten[%]/.min[[2]]]

min2={b->b2}/.%//N

nint2= NIntegrate[(Subscript[Ψ, 3])[x]^2/.a->1/.min2,{x,0,1}]

Plot[{Sqrt[2]Sin[3Pi x],(Subscript[Ψ, 3])[x]/Sqrt[nint]/.{a->1,b->-4.633}},{x,0,1}]

Eb[b]/Pi^2/.min2

Our trial wavefunction for the second excited state isn' t as good as the one for the ground state.  It hasn' t done badly on the energy though (10.3 instead of 9).  In this case it is greater than the true value because it is so nearly orthogonal to the true ground state, but this wasn't guaranteed.

We can look at the expansion of this state in terms of the true states of the system (since here we do actually know them) - we will set a-> 1  for numerical integration:

coef=Table[Sqrt[2/nint2] NIntegrate[Sin[n Pi x](Subscript[Ψ, 3])[x]/.a->1/.min2,{x,0,1}],{n,1,41,2}]

So the overlap with the ground state really is tiny.  We can also check that <H>=(∑_n)(Subscript[P, n])(Subscript[E, n])

Sum[coef[[n]]^2,{n,1,21}]

Sum[coef[[n]]^2( 2n-1)^2,{n,1,21}]

With 21 non-vanishing terms, we have nearly converged on the norm and on  <H>: