Computing the expectation value of \(\langle x^4 \rangle\) for the 1-dimensional particle-in-a-box

Here is an example of calculating \(\langle x \rangle\) and \(\langle x^2 \rangle\) for the 1-dimensional particle-in-a-box.

The average position is expected to be

\[ \langle x \rangle =\tfrac{a}{2} \]

We can confirm this by explicit integration,

\[\begin{split} \begin{align} \langle x \rangle &= \int_0^a \psi_n^*(x)\, x \,\psi_n(x) dx \\ &= \int_0^a \left(\sqrt{\tfrac{2}{a}} \sin\left(\tfrac{n \pi x}{a} \right)\right) x \left(\sqrt{\tfrac{2}{a}}\sin\left(\tfrac{n \pi x}{a} \right)\right) dx \\ &= \tfrac{2}{a} \int_0^a x \sin^2\left(\tfrac{n \pi x}{a} \right) dx \\ &= \tfrac{2}{a} \left[ \tfrac{x^2}{4} - x \tfrac{\sin \tfrac{2n \pi x}{a}}{\tfrac{4 n \pi}{a}} - \tfrac{\cos \tfrac{2n \pi x}{a}}{\tfrac{8 n^2 \pi^2}{a^2}} \right]_0^a \\ &= \tfrac{2}{a} \left[ \tfrac{a^2}{4} - 0 - 0 \right] \\ &= \tfrac{a}{2} \end{align} \end{split}\]

Similarly, we expect that the expectation value of \(\langle x^2 \rangle\) will be proportional to \(a^2\). We can confirm this by explicit integration,

\[\begin{split} \begin{align} \langle x^2 \rangle &= \int_0^a \psi_n^*(x)\, x^2 \,\psi_n(x) dx \\ &= \int_0^a \left(\sqrt{\tfrac{2}{a}} \sin\left(\tfrac{n \pi x}{a} \right)\right) x^2 \left(\sqrt{\tfrac{2}{a}}\sin\left(\tfrac{n \pi x}{a} \right)\right) dx \\ &= \tfrac{2}{a} \int_0^a x^2 \sin^2\left(\tfrac{n \pi x}{a} \right) dx \\ &= \tfrac{2}{a} \left[ \tfrac{x^3}{6} - x^2 \tfrac{\sin \tfrac{2n \pi x}{a}}{\tfrac{4 n \pi}{a}} - x \tfrac{\cos \tfrac{2n \pi x}{a}}{\tfrac{4 n^2 \pi^2}{a^2}} - \tfrac{\sin \tfrac{2n \pi x}{a}}{\tfrac{8 n^3 \pi^3}{a^3}} \right]_0^a \\ &= \tfrac{2}{a} \left[ \tfrac{a^3}{6} - 0 - \tfrac{a}{{\tfrac{4 n^2 \pi^2}{a^2}}} - 0 \right] \\ &= \tfrac{2}{a} \left[ \tfrac{a^3}{6} - \tfrac{a^3}{4 n^2 \pi^2} \right] \\ &= a^2\left[ \tfrac{1}{3} - \tfrac{1}{2 n^2 \pi^2} \right] \end{align} \end{split}\]

Task:

Use the above example to calculate \(\langle x^4 \rangle\). You can follow next steps:

1. Complete the compute_wavefunction function by coding the expression for calculating the wavefunction. Don’t forget the normalization constant. Also, take into account that \(x\) values could be negative and greater than \(a\), and that \(x\) values could be a simple real number (a float) or a list of numbers (a numpy.ndarray).

2. Using the compute_wavefunction function, complete the compute_probability function.

3. Fill the gaps in compute_moments function; this is especially easy if you use the compute_probability function.

4. Complete the function check_moment.

5. Now you are ready to calculate the average. Use the quad function for numerical integration and complete the calc_average function.

6. Use your functions to calculate the \(\langle x^2 \rangle\) and compare it with the analytic value. We encourage you to play with different values of \(n\) and \(a\).

7. Submit your notebook as described at the first tutorial.

import numpy as np
from scipy.integrate import quad



# Define a function for the wavefunction
def compute_wavefunction(x, n, a):
    """Compute 1-dimensional particle-in-a-box wave-function value(s).

    Parameters
    ----------
    x: float or np.ndarray
        Position of the particle.
    n: int
        Quantum number value.
    a: float
        Length of the box.
    """
    # check argument n
    if not (isinstance(n, int) and n > 0):
        raise ValueError("Argument n should be a positive integer.")
    # check argument a
    if a <= 0.0:
        raise ValueError("Argument a should be positive.")
    # check argument x
    if not (isinstance(x, float) or hasattr(x, "__iter__")):
        raise ValueError("Argument x should be a float or an array!")

    # compute wave-function

    ### START YOUR CODE HERE
    value = np.sqrt(2/a)*np.sin(np.pi*x*n/a)
    ### END YOUR CODE HERE


    # set wave-function values out of the box equal to zero

    ### START YOUR CODE HERE
    if hasattr(x, "__iter__"):
      value[np.logical_or(x>a, x<0)] = 0
    elif x<0 or x>a:
      value = 0

    ### END YOUR CODE HERE
    return value



# Define a function for the wavefunction squared
def compute_probability(x, n, a):
    """Compute 1-dimensional particle-in-a-box probablity value(s).

    See `compute_wavefunction` parameters.
    """
    ### START YOUR CODE HERE
    probability = compute_wavefunction(x, n, a)**2
    ### END YOUR CODE HERE
    return probability



def compute_moment(x, n, a, power):
    """Compute the x^power moment of the 1-dimensional particle-in-a-box.

    See `compute_wavefunction` parameters.
    """
    return compute_probability(x, n, a)*x**power



#Compute <x^power>, the expectation value of x^power
def calc_average(n, a, power):
    """
    Compute the average value by explicit integrating
    """

    ### START YOUR CODE HERE
    avg, error = quad(compute_moment, 0, a, args=(n, a, power))
    ### END YOUR CODE HERE

    return avg


#This next bit of code just prints out the values.
def check_moments(n, a):
    #check the computed values of the moments against the analytic formula

    ### START YOUR CODE HERE
    power = 2
    avg_r2 = calc_average(n, a, power)
    ### END YOUR CODE HERE

    print(f"<r^2> computed = {avg_r2:.5f}")
    print(f"<r^2> analytic = {a**2*(1/3 - 1./(2*n**2*np.pi**2))}")


#Principle quantum number:
n = 1

#Box length:
a = 1

check_moments(a, n)

<r^2> computed = 0.28267
<r^2> analytic = 0.2826727415121644