Variation Method for the Particle in a Box.
Particle in an infinite square well.
Expand in a basis of polynomials in x.


Statement of the Problem

A particle with mass meff is confined to move in one dimension under the potential V(x) = {0 if -L/2 < x < L/2, otherwise}. Schroedinger's equation for this system is
wpe13.gif (2141 bytes)

Use the variation principle to find approximate eigenvalues and eigenfunctions for a trial function having the form of a polynomial summation.

Dimensionless Form

With a little effort the Schroedinger equation can be transformed into an equivalent but simpler equation without units. The new dimensionless Schroedinger equation applies to all PIB systems regardless of the particle mass or the length of the box. First, remove the units of length by defining x=yL/2 for the potential is zero when -1<y<1 and infinte outside this range. Thus, the box is referenced to dimensionless coordinate y and has been standardized to length 2 with the center at y=0. This change of coordinates induces a change of the wave function: u(y) = y(yL/2). Second, introduce the "natural" or characteristic energy wpe16.gif (1224 bytes)   and divide both sides of Schroedinger's equation by this amount. The result is the following dimensionless Schroedinger equation:
wpe17.gif (2071 bytes),
where e=E/Enatural.

Notice that the dimensionless Schroedinger equation contains no memory of the particle mass or the box length - this information is contained in the coordinate transformation x=yL/2 and the natural energy. The goal will be to solve the dimensionless form of the problem and then restore units for the particular mass and box at hand.

Choosing the Basis Functions

It is crucial to satisfy the boundary conditions: any basis functions used must satisfy the boundary conditions. With this in mind, some possible variational trial functions that might be part of a basis are f1(y)=(y-1)(y+1), f2(y)=cos(py/2), and f3(y)=cosh(y)-cosh(1). In order to perform the linear variation we require a large set of linearly independent trial functions.


As an example, here is such a set generalized from f1:
wpe19.gif (1306 bytes), with n=0,1,...,M.
For instance, when M=2 the three basis functions are f0=(y-1)3(y+1), f1=(y-1)2(y+1)2, and f2=(y-1)(y+1)3. Ultimately, the basis functions are chosen for reasons of simplicity, convenience, and aesthetics.

Basis functions are only useful if the required overlap and hamiltonian matrix elements (integrals) can be readily computed. We illustrate using the example basis.

Overlap and Hamiltonian Matrices

Elements of the overlap matrix are computed from the definition wpe1A.gif (1286 bytes).

For the polynomial basis functions of the example, this is straightforward though tedious. Both basis functions can be expressed in a binomial expansion the the powers of y can be integrated term by term. From the definition, it follows that the overlap matrix has to be symmetric: Si,j = Sj,i for our example.

Elements of the hamiltonian matrix are computed from wpe1B.gif (1608 bytes). This time, the second derivative of basis functions has to be computed and again integrals performed. For our example, the derivatives are polynomials in y so that the matrix element evaluation is straightforward. As a check, the hamiltonian matrix also is symmetric: Hi,j = Hj,i.

Numerical values of S and H are shown in the complete solution file.

Secular equation and secular determinant

Eigenvalues are found from the secular determinant: |H-WS|=0, a M-th order polynomial with M roots. Denote these roots by W0<W1<...WM.   For the ground state, the smallest root, W0, is used. Higher roots are bounds for excited states.

For the example, with M=3 (four basis functions) the lowest root is 1.00001471 (an error of only ~15 ppm from the exact ground state energy).

Wave functions are found by solving the secular equation: (H-WiS)ai=0, and substituting the eigenvector into the basis expansion: wpe1C.gif (1282 bytes). Note that the vector ai is labelled for its eigenvalue and its elements also carry a label for the basis function. The collection of ai,j comprise a square matrix the columns of which are the the eigenvectors belonging to Wi.

For the example, see the calculated S and H matrices in the solution file.

Error analysis

Errors in the variational approximation can be assessed in several ways.

These methods are illustrated for our example.

Example illustrations:

1. compare trial energy with exact energy and
assess convergence of energy with basis size.

Here is a table to summarize the calculated variational energy values as the basis set is expanded. Recall there are 2M+1 basis functions. The table illustrates a general feature of linear variation: Convergence is usually fastest for the ground state and becomes slower with excitation (presumeably it becomes harder to fit the nodes in multiply excited wave function).

M (2M+1)  1 (3) 3  (7) 5 (11) 7 (15)
W0   1.013211836 1.000014714   1.000000011   1.000000004
W1   4.255489713 4.002344086   4.000005654   4.000000013
W2   - 10.347957854   9.035177135   9.000306253
W3   - 20.314739988   16.210526743   16.004427334
W4   - -   35.559363311   25.777803326
W5    - - 57.805391048   38.148049011


- -   - 89.048113554
W7   - -   -   131.556347522

2. Wave functon comparison.

Ground state wave functions are compared in the following graphs. The first graph shows the approximate variational wave function (F(x)) for M=3 and the exact wave function. To the eye, there is no difference between approximate and exact solutions. The second graph shows the relative error of the trial function. (The sign of the trial function was reversed in these plots. Overall wave function sign (phase) has no physical significance.) Here the relative error of the trial function are obvious; it is largest near the ends of the box where y(y) is small and less than about 2 ppt throughout most of the box.

wpe20.gif (5174 bytes)

3. Local energy function

According to the Schroedinger equation, the quotient function defined as Hy(x)/y(x) = E is a constant. But for a trial function the corresponding quotient, e(x) = Hf(x)/f(x) may vary across the box. This function is called the "local energy" function. Any departure from constant value reveals an error in the trial function. For example, when M=3 the local energy function is graphed below.

wpe1F.gif (3674 bytes)

For the dimensionless problem this local energy deviates from Etrue(=1) most greatly near the ends of the box and stays within about 3% of Etrue through most of the box. [This solution was created using mathcad. The solution file is available for downloading. It can easily be extended to a larger basis. Minor changes in the mathcad solution allow the potential to be "decorated"; e.g., a central well (or barrier) can be added to the box.]

4. Lower bounds

There are other ways to assess the errors in a trial variational function. Since the variation principle provides an upper bound to the true energy, it is natural to seek a lower bound to bracket the exact energy. Such lower bounds can be found but often they exagerate the size of the error. One such lower bound method is due to Temple [G. Temple, Proc. Roy. Soc., A119, 276 (1928)] and involves computing matrix elements of the square of the hamiltonian (a computationally laborious process).


  1. For the example problem, plot the trial wave function for the first excited state in the box, and the relative error of this state function.
  2. Trial functions don't have to be polynomials. Calculate the variational energy from the trial function f3(y)=cosh(ay)-cosh(a) wherea  is a variable parameter. Hint: The Rayleigh quotient is indeterminate at a=0. Make a graph of W(a) versus a (near but not including a=0) to estimate the minimum of W. Use L'Hospital's rule to determine the minimum. How big is the error in this trial energy?
  3. Symmetry was ignored in the treatment of the example on this page. To consider symmetry, answer the following:
    (a) Show that fn(-y) = (-1)M fM-n(y) . Hence, also that fM-n(-y)= (-1)Mfn(y). Finally, show that the sum and differences fn(y) fM-n(y) are either even parity or odd parity depending on M.
    (b) Examine the linear coefficients of the ground state trial function in the example. Namely, consider the expansion coefficients a0 = -0.166284, a1 = 0.687277, a2 = -0.687277, a3 = 0.166284 for M=3. Explain the magnitudes and signs in light of symmetry.
  4. Use the following basis set for a linear variational treatment of the PIB: fn(y) = (y2-1)n, n=1,2,...,M. (a) Calculate the variational energy and wave function of the ground state for M=1,3,5. (b) Analyse the errors in these variational results.


Created or up-dated 08/03/99   by R.D. Poshusta
HH01580A.gif (1311 bytes)