# Computational Exercise 1: A Vary Great, Very Harmonic, Wavefunction
We will apply the Variational Method to the Harmonic Oscillator.  In particular, we will optimize
the trial wavefunction given by 
\begin{equation}
\Phi(x) = {\rm exp}\left( -\alpha x^2 \right),
\end{equation}
where $\alpha$ is an adjustable parameter that will be variationally optimized.  In particular, the goal
will be to minimize the energy functional given by
\begin{equation}
E[\Phi(x)] = \frac{\int_{-\infty}^{\infty} \Phi^* (x) \: \hat{H} \: \Phi(x) dx }{\int_{-\infty}^{\infty} \Phi^* (x) \: \Phi(x) dx }
\end{equation}
by modulating the trial wavefunction through the adjustable parameter $\alpha$.

The Hamiltonian operator for Harmonic Oscillator is given by 
\begin{equation}
\hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \frac{1}{2} k \: x^2.
\end{equation}
We will use a system of natural units where $\hbar = 1$, the mass of the particle $m = 1$, and the force constant of
the potential $k = 1$.

### Questions
1. What is the optimal value of $E[\Phi(x)]$ in natural units?  How does this value compare to the true ground state energy of the Harmonic Oscillator?
2. What is the optimal value of $\alpha$ in natural units?  How does $\Phi(x)$ with this value of $\alpha$ compare to the true ground state wavefunction?
3. What is the kinetic energy expectation value of the optimized trial wavefunction in natural units?
4. What is the potential energy expectation value of the optimized trial wavefunction in natural units?

### (Hopefully) Helpful Hints
You can represent your trial wavefunction numerically on the computer by creating an array of x-values, evaluating
the formula for your trial wavefunction at each of those x-values, and then storing each of those function values in another array.  The following lines of code will do this using $numpy$ arrays and the $numpy.exp$ function. 



```python
import numpy as np
from matplotlib import pyplot as plt
### create an array of 1000 x-values between x=-20 and x=20 atomic units
x_array = np.linspace(-20,20,1000)
### pick a trial value of alpha... might as well be 1 atomic unit!
alpha = 1.0
### evaluate your trial wavefunction at each of the x-values stored in the array called 'x_array'
Phi = np.exp(-alpha * x_array**2)
Vx = 1/2. * x_array **2

plt.plot(x_array, Phi, 'red', label='Trial')
plt.plot(x_array, Vx, 'green', label='Potential')
plt.xlim(-5,5)
plt.ylim(0,10)
plt.legend()
plt.show()


```

To evaluate the kinetic energy of a trial wavefunction, we will need to take the second derivative of the trial wavefunction.
One way to take derivatives on a computer is to use the method of finite differences.  In particular, the forward finite difference
approximation to the first derivative of your trial wavefunction $\Phi(x)$ at the point $x_0$ can be written symbolically as
\begin{equation}
\frac{d \Phi(x_0)}{dx} = \frac{\Phi(x_0 + dx) - \Phi(x_0)}{dx}
\end{equation}
where in our code, $dx$ will be determined by the range of x-values and the number of x-values used in defining
$x_{array}$.

Create a function called $dfdx$ that takes an array of x-values and an array of function values evaluated at each x-value and
returns the first derivative of the function values with respect to x.  Then, show that calling this function in sequence can be used to compute the second derivative of an array of function values.  

A skeleton of the function looks like this:


```python
def dfdx(fx, x):
    ### loop through values of fx and x and use forward finite-differences 
    ### to compute the d/dx of fx at each point and store in an array called 'fxp'
    ### make sure to poperly handle
    ### the last point where f(x_0 + dx) is not stored
    ### create a numpy array called fxp that has the same length
    ### and data type as the array of function values you are differentiating
    fxp = np.zeros_like(fx)
    
    for i in range(0,len(x)-1):
        rise = fx[i+1]-fx[i]
        run = x[i+1]-x[i]
        slope = rise/run
        fxp[i] = slope

    ### once you have computed all fxp values, return the array!
    return fxp
```

And the calling sequence for $dfdx$ to compute the second derivative of a trial wavefunction we defined previously will look like this:


```python
### Call dfdx on Phi to get the first derivative of Phi
Phi_p = dfdx(Phi, x_array)

Phi_wa = -2 * alpha * x_array * Phi

plt.plot(x_array, Phi_p, 'red', label='dfdx result')
plt.plot(x_array, Phi_wa, 'b--', label='Wolfram Alpha result')
plt.legend()
plt.show()


### Call dfdx on Phi_p to get the second derivative of Phi
Phi_pp = dfdx(Phi_p, x_array)
```


```python

```

ae70e694f18c74fe912c32b074c4af03d43f858f
ipynb
Jupyter Notebook
notebooks/Harmonic Oscillator.ipynb
FoleyLab/FoleyLab.github.io
1f84e4dc2f87286dbd4e07e483ac1e48943cb493
CC-BY-3.0
notebooks/Harmonic Oscillator.ipynb
FoleyLab/FoleyLab.github.io
1f84e4dc2f87286dbd4e07e483ac1e48943cb493
CC-BY-3.0
2020-02-25T08:45:19.000Z
2021-05-19T04:28:30.000Z
notebooks/Harmonic Oscillator.ipynb
FoleyLab/FoleyLab.github.io
1f84e4dc2f87286dbd4e07e483ac1e48943cb493
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES