The variational principle and the quantum mechanical particle-in-a-box

Consider a one-dimensional quantum mechanical particle in a box ( $-1\le x \le 1$ ) described by the Schrödinger equation:

-\frac{\mathrm{d}^2\psi}{\mathrm{d}x^2} = E\psi,

in energy units for which $\hbar^2/(2m) = 1$ with $$m$$ the mass of the particle. The exact solution for the ground state of this system is given by

\psi = \cos\left(\frac{\pi x}{2}\right), \quad E = \frac{\pi^2}{4}.

An approximate solution may be arrived at using the variational principle by minimizing the expectation value of the energy of a trial wavefunction,

\psi_\mathrm{trial} = \sum_{n=0}^N a_n\phi_n(x)

with respect to the coefficients $$a_n$$ . Taking the basis functions to have the following symmetrized polynomial form,

\phi_n = (1-x)^{N-n+1}(x+1)^{n+1},

use scipy.optimize.minimize and scipy.integrate.quad to find the optimum value of the expectation value (Rayleigh-Ritz ratio):

\mathcal{E} = \frac{\langle \psi_\mathrm{trial} | \hat{H} | \psi_\mathrm{trial} \rangle}{\langle \psi_\mathrm{trial} | \psi_\mathrm{trial} \rangle} = -\frac{\int_{-1}^1\psi_\mathrm{trial}\frac{\mathrm{d}^2}{\mathrm{d}x^2}\psi_\mathrm{trial}\;\mathrm{d}x}{\int_{-1}^1\psi_\mathrm{trial}\psi_\mathrm{trial}\;\mathrm{d}x}.

Compare the estimated energy, $\mathcal{E}$ , with the exact answer for $$N=1,2,3,4$$ .

[Hint: Use np.polynomial.Polynomial objects to represent the basis and trial wavefunctions].

The code below minimizes the expectation value of the energy for the polynomial approximation to the particle-in-a-box wavefunction. The list basis holds the basis set as Polynomial objects, in which the trial wavefunction is expanded with coefficients held in the array a. The necessary integrals are carried out by functions S and H using scipy.integrate.quad.

variationalpib

import sys
import numpy as np
from scipy.integrate import quad
from scipy.optimize import minimize

Polynomial = np.polynomial.Polynomial
import matplotlib.pyplot as plt


def get_basis(N):
    """Return the N+1 basis functions for n = 0, 1, ..., N."""
    p1, p2 = Polynomial([1, -1]), Polynomial([1, 1])
    basis = [p1 ** (N - n + 1) * p2 ** (n + 1) for n in range(N + 1)]
    return basis


def get_psi(a, basis):
    """
    Return a Polynomial object representing the trial wavefunction expanded
    in basis with coefficients a.

    """
    return sum([a[i] * p for i, p in enumerate(basis)])


def S(psi):
    """Return the value of the normalization integral, <psi|psi>."""
    norm, _ = quad(psi**2, 1, -1)
    return norm


def H(psi):
    """Return the value of the integral <psi | H | psi>."""
    d2psi = psi.deriv(2)
    integral, _ = quad(psi * d2psi, 1, -1)
    return integral


def rayleigh_ritz(a, basis):
    """Return the value of the Rayleigh-Ritz ratio."""
    psi = get_psi(a, basis)
    return -H(psi) / S(psi)


fig, axes = plt.subplots(nrows=2, ncols=2)
x = np.linspace(-1, 1, 1000)

for N in range(1, 5):
    print("-" * 60, f"N = {N}", sep="\n")
    basis = get_basis(N)
    # Initial guess for the coefficients in the trial wavefunction
    a0 = np.ones(N + 1)

    # Minimize the Rayleigh-Ritz ratio w.r.t. a
    res = minimize(rayleigh_ritz, a0, args=(basis,))
    if not res["success"]:
        print("Variational minimization failed:")
        print(res)
        sys.exit(1)

    # Report and compare with the exact answer.
    print("Variational optimum E =", res["fun"])
    print("     Exact energy, E0 =", np.pi**2 / 4)

    # Plot the true and approximate wavefunctions
    psi = get_psi(res["x"], basis)
    # Obtain the approximate wavefunction and normalize it, ensuring it's
    # positive.
    psi_approx = abs(psi(x)) / np.sqrt(abs(S(psi)))
    # Exact ground-state wavefunction
    psi_exact = np.cos(np.pi * x / 2)

    ax = axes[(N - 1) // 2, (N - 1) % 2]
    ax.plot(x, psi_approx, label="Approximation", lw=1)
    ax.plot(x, psi_exact, label="Exact", lw=1)
    ax.legend(loc="lower center", fontsize=10)
    ax.set_title(rf"$N={N:1d}$")

plt.tight_layout()
plt.show()

We also plot the exact and approximate wavefunctions, but they are indistinguishable for $$N>1$$ at this scale. The output shows good agreement of the exact and estimated energies for $$N=4$$ :

------------------------------------------------------------
N = 1
Variational optimum E = 2.500000000000001
     Exact energy, E0 = 2.4674011002723395
------------------------------------------------------------
N = 2
Variational optimum E = 2.467437405333598
     Exact energy, E0 = 2.4674011002723395
------------------------------------------------------------
N = 3
Variational optimum E = 2.467437405329251
     Exact energy, E0 = 2.4674011002723395
------------------------------------------------------------
N = 4
Variational optimum E = 2.4674011087466323
     Exact energy, E0 = 2.4674011002723395

Note that the variational optimum is greater than the true ground state energy, as required by the variational principle.

Approximate variational solutions to the particle-in-a-box problem with a polynomial basis.

P8.4.4The variational principle and the quantum mechanical particle-in-a-box