P8.4.4: The variational principle and the quantum mechanical particle-in-a-box
Question P8.4.4
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].
Solution P8.4.4
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.
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.
