# 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{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} 
$$

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{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} 
$$

### 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](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) 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](Tutorial_0.ipynb).


```python
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 YOU 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 YOU CODE HERE
    if hasattr(x, "__iter__"):
      value[np.logigal_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


d66cacda4a3a2e8f5149ae6ca3194942bc0e16e5
ipynb
Jupyter Notebook
book/x4_mocked.ipynb
guacke/IntroQM2022
3a268eda86491db86f674784169e921ff6a22b2b
CC0-1.0
2022-02-08T18:42:37.000Z
2022-02-21T19:33:46.000Z
book/x4_mocked.ipynb
guacke/IntroQM2022
3a268eda86491db86f674784169e921ff6a22b2b
CC0-1.0
2022-01-26T18:45:29.000Z
2022-03-04T20:32:52.000Z
book/x4_mocked.ipynb
guacke/IntroQM2022
3a268eda86491db86f674784169e921ff6a22b2b
CC0-1.0
2022-02-08T17:56:55.000Z
2022-03-03T08:30:53.000Z
Qwen/Qwen-72B
1. YES
2. YES