# The Harmonic Oscillator Strikes Back

*Note:* Much of this is adapted/copied from https://flothesof.github.io/harmonic-oscillator-three-methods-solution.html

This week we continue our adventures with the harmonic oscillator.

 

The harmonic oscillator is a system that, when displaced from its equilibrium position, experiences a restoring force F proportional to the displacement x:

$$F=-kx$$

The potential energy of this system is 

$$V = {1 \over 2}k{x^2}$$

These are sometime rewritten as

$$ F=- \omega_0^2 m x, \text{ }  V(x) = {1 \over 2} m \omega_0^2 {x^2}$$

Where $\omega_0  = \sqrt {{k \over m}} $

If the equilibrium value of the harmonic oscillator is not zero, then 

$$ F=- \omega_0^2 m (x-x_{eq}), \text{ }  V(x) = {1 \over 2} m \omega_0^2 (x-x_{eq})^2$$

## 1. Harmonic oscillator from last time (with some better defined conditions)

Applying the harmonic oscillator force to Newton's second law leads to the following second order differential equation

$$ F = m a $$

$$ F= -m \omega_0^2 (x-x_{eq}) $$

$$ a = - \omega_0^2 (x-x_{eq}) $$

$$ x(t)'' = - \omega_0^2 (x-x_{eq}) $$

The final expression can be rearranged into a second order homogenous differential equation, and can be solved using the methods we used above

This is already solved to remind you how we found these values


```python
import sympy as sym
sym.init_printing()
```

**Note** that this time we define some of the properties of the symbols. Namely, that the frequency is always positive and real and that the positions are always real


```python
omega0,t=sym.symbols("omega_0,t",positive=True,nonnegative=True,real=True)
xeq=sym.symbols("x_{eq}",real=True)
x=sym.Function("x",real=True)
x(t),omega0
```


```python
dfeq=sym.Derivative(x(t),t,2)+omega0**2*(x(t)-xeq)
dfeq
```


```python
sol = sym.dsolve(dfeq)
sol
```


```python
sol,sol.args[0],sol.args[1]
```

**Note** this time we define the initial positions and velocities as real


```python
x0,v0=sym.symbols("x_0,v_0",real=True)
ics=[sym.Eq(sol.args[1].subs(t, 0), x0),
     sym.Eq(sol.args[1].diff(t).subs(t, 0), v0)]
ics
```


```python
solved_ics=sym.solve(ics)
solved_ics
```

### 1.1 Equation of motion for $x(t)$


```python
full_sol = sol.subs(solved_ics[0])
full_sol
```

### 1.2 Equation of motion for $p(t)$


```python
m=sym.symbols("m",positive=True,nonnegative=True,real=True)
p=sym.Function("p")
sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t))
```

## 2. Time average values for a harmonic oscillator

If we want to understand the average value of a time dependent observable, we need to solve the following integral


$${\left\langle {A(t)} \right\rangle}_t  = \begin{array}{*{20}{c}}
{\lim }\\
{\tau  \to 0}
\end{array}\frac{1}{\tau }\int\limits_0^\tau  {A(t)dt} $$

### 2.1 Average position ${\left\langle {x} \right\rangle}_t$ for a harmonic oscillator


```python
tau=sym.symbols("tau",nonnegative=True,real=True)
xfunc=full_sol.args[1]
xavet=(xfunc.integrate((t,0,tau))/tau).limit(tau,sym.oo)
xavet
```

The computer does not always make the best choices the first time. If you treat each sum individually this is not a hard limit to do by hand. The computer is not smart. We can help it by inseting an `expand()` function in the statement


```python
xavet=(xfunc.integrate((t,0,tau))/tau).expand().limit(tau,sym.oo)
xavet
```

### 2.2 Excercise: Calculate the average momenta  ${\left\langle {p} \right\rangle}_t$ for a harmonic oscillator


```python
import sympy as sym
sym.init_printing()
```


```python
m=sym.symbols("m",positive=True,nonnegative=True,real=True)
p=sym.Function("p")
sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t))
```


```python
tau=sym.symbols("tau",nonnegative=True,real=True)
pfunc=sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t)).args[1]
pavet=(pfunc.integrate((t,0,tau))/tau).limit(tau,sym.oo)
pavet
```

### 2.3 Exercise: Calculate the average kinetic energy of a harmonic oscillator


```python
kefunc=(sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t)).args[1])**2/(2*m)
keavet=(kefunc.integrate((t,0,tau))/tau).expand().limit(tau,sym.oo)
keavet
```


```python
sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t)).args[1]**2/(2*m)
```


```python
((m*sol.args[1].subs(solved_ics[0]).diff(t)**2/(2*m)))
```


```python
m*sol.args[1].subs(solved_ics[0]).diff(t)**2/(2*m)
```

## 3.  Ensemble (Thermodynamic) Average values for a harmonic oscillator

If we want to understand the thermodynamics ensemble average value of an observable, we need to solve the following integral. 


$${\left\langle {A(t)} \right\rangle}_{T}  = \frac{\int{A e^{-\beta H}dqdp}}{\int{e^{-\beta H}dqdp} } $$

You can think of this as a Temperature average instead of a time average.

Here $\beta=\frac{1}{k_B T}$ and the classical Hamiltonian, $H$ is

$$ H = \frac{p^2}{2 m} + V(q)$$

**Note** that the factors of $1/h$ found in the classical partition function cancel out when calculating average values

### 3.1 Average position ${\left\langle {x} \right\rangle}_t$ for a harmonic oscillator

For a harmonic oscillator with equilibrium value $x_{eq}$, the Hamiltonian is 
$$ H = \frac{p^2}{2 m} + \frac{1}{2} m \omega_0 (x-x_{eq})^2 $$

First we will calculate the partition function $\int{e^{-\beta H}dqdp}$


```python
k,T=sym.symbols("k,T",positive=True,nonnegative=True,real=True)
xT,pT=sym.symbols("x_T,p_T",real=True)
ham=sym.Rational(1,2)*(pT)**2/m + sym.Rational(1,2)*m*omega0**2*(xT-xeq)**2
beta=1/(k*T)
bolz=sym.exp(-beta*ham)
z=sym.integrate(bolz,(xT,-sym.oo,sym.oo),(pT,-sym.oo,sym.oo))
z
```

Then we can calculate the numerator $\int{A e^{-\beta H}dqdp}$



```python
numx=sym.integrate(xT*bolz,(xT,-sym.oo,sym.oo),(pT,-sym.oo,sym.oo))
numx
```

And now the average value


```python
xaveT=numx/z
xaveT
```

### 3.2 Exercise: Calculate the average momenta  ${\left\langle {p} \right\rangle}_t$ for a harmonic oscillator

After calculating the value, explain why you think you got this number 


```python
nump=sym.integrate(pT*bolz,(xT,-sym.oo,sym.oo),(pT,-sym.oo,sym.oo))
nump
```


```python
paveT=nump/z
paveT
```

### 3.3 Exercise: Calculate the average kinetic energy

The answer you get here is a well known result related to the energy equipartition theorem


```python
keaveT=sym.integrate(pT**2/(2*m)*bolz,(xT,-sym.oo,sym.oo),(pT,-sym.oo,sym.oo))/z
keaveT
```

# Back to the lecture

## 4. Exercise Verlet integrators

In this exercise we will write a routine to solve for the equations of motion for a hamonic oscillator.

Plot the positions and momenta (seprate plots) of the harmonic oscillator as a functions of time.

Calculaate trajectories using the following methods:
1. Exact solution
2. Simple taylor series expansion
3. Predictor-corrector method
4. Verlet algorithm
5. Leapfrog algorithm
6. Velocity Verlet algorithm


```python
# Your code here
```


```python
%matplotlib inline
```


```python
omega0,t=sym.symbols("omega_0,t",positive=True,nonnegative=True,real=True)
xeq=sym.symbols("x_{eq}",real=True)
x=sym.Function("x",real=True)
full_sol1 = sym.simplify(full_sol.subs({x0:10, xeq:0 , v0:10, omega0:1}))
sym.plot(full_sol1.rhs,(t,-10,10))
```


```python
m=sym.symbols("m",positive=True,nonnegative=True,real=True)
p=sym.Function("p")
sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t))
momentum=sym.Eq(p(t),m*sol.args[1].subs(solved_ics[0]).diff(t))
```


```python
momentum1=sym.simplify(momentum.subs({x0:10, xeq:0, v0:10, omega0:1, m:1}))
```


```python
sym.plot(momentum1.rhs,(t,-10,10))
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=0
t=0.25
v_t0=1
xlist=[]
for i in range(0,500):
    v_t=v_t0+(1/2)*t**2*(-(x_t0))
    v_t0=v_t
    x_t= x_t0+(v_t0)*(t)
    x_t0=x_t
    xlist.append(x_t)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Taylor series approximation')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=0
x_t1=1
t=0.4
a=1
xlist=[]
for i in range(0,100):
    x_t2=2*x_t1-x_t0+t**2*(-(x_t1))
    x_t0=x_t1
    x_t1=x_t2
    xlist.append(x_t2)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Verlet Algorithm-Position')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=0
x_t2=1
t=2
xlist=[]
for i in range(0,100):
    v_t=(x_t2-x_t0)/2*t
    x_t0=x_t2
    x_t2=v_t
    xlist.append(v_t)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Verlet Algorithm-Velocity')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=1
v_minushalft=0
t=0.2
xlist=[]
for i in range(0,100):
    v_halft=v_minushalft+(t)*(-(x_t0))
    v_minushalft=v_halft
    x_t1=x_t0+(t)*(v_halft)
    x_t0=x_t1
    xlist.append(x_t1)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Leapfrog Algorithm-Position')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=1
t=0.3
v_minushalft=1
v_halft=2
xlist=[]
for i in range(0,100):
    v_t=(1/2)*((v_halft)+(v_minushalft))
    v_minushalft=v_t
    v_halft=v_minushalft+(t)*(-(x_t0))
    v_minushalft=v_halft
    x_t1=x_t0+(t)*(v_halft)
    x_t0=x_t1
    xlist.append(v_t)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Leapfrog Algorithm-Velocity')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
dt=0.1
x_0=1
v_0=1
xlist=[]
for i in range(0,100):
    x_1 = x_0 + v_0 * dt + 0.5 * -x_0 * dt * dt
    x_0 = x_1
    v_1 = v_0 + 0.5 * (-x_0 + -x_1) * dt
    v_0 = v_1
    xlist.append(x_1)
    #print(xlist)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Velocity Verlet-Position Approx')
```


```python
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
x_t0=1
t=2
v_t0=1
xlist=[]
for i in range(0,50):
    v_halft=v_t0+(1/2)*t*(-x_t0)
    x_t0=v_halft
    xlist.append(v_halft)
plt.plot(xlist)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.title('Velocity Verlet-Velocity Approximation')
```


```python
#i officially quit and resign
```


```python

```

bfcf4bf1d92094663e973086615f40eeaa46149e
ipynb
Jupyter Notebook
harmonic_student.ipynb
sju-chem264-2019/10-14-19-more-harmonic-oscillator-deannapatti
1632d729fd7277cd9a2bee4a4f1b8217294a9e12
MIT
harmonic_student.ipynb
sju-chem264-2019/10-14-19-more-harmonic-oscillator-deannapatti
1632d729fd7277cd9a2bee4a4f1b8217294a9e12
MIT
harmonic_student.ipynb
sju-chem264-2019/10-14-19-more-harmonic-oscillator-deannapatti
1632d729fd7277cd9a2bee4a4f1b8217294a9e12
MIT
Qwen/Qwen-72B
1. YES
2. YES