# 5.24. Computing the expectation value of $$\langle x^4 \rangle$$ for the 1-dimensional particle-in-a-box¶

## 5.24.1. 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}

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

# 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 YOU CODE HERE
value = np.sqrt(2/a)*np.sin(np.pi*x*n/a)

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

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

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.
"""
probability = compute_wavefunction(x, n, a)**2
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
"""

avg, error = quad(compute_moment, 0, a, args=(n, a, power))

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

power = 2
avg_r2 = calc_average(n, a, power)

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