# 10 ODE integrators: Verlet


```python
from importlib import reload

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.style.use('ggplot')

import integrators2_solution as integrators2
```


```python
reload(integrators2)
```




    <module 'integrators2_solution' from '/Volumes/Data/oliver/Documents/Teaching/ASU/CompPhys_PHY494/2018/PHY494-resources/10_ODEs/integrators2_solution.py'>



## Velocity Verlet

Use expansion *forward* and *backward* (!) in time (Hamiltons (i.e. Newton without friction) equations are time symmetric)

\begin{align}
r(t + \Delta r) &\approx r(t) + \Delta t\, v(t) + \frac{1}{2m} \Delta t^2 F(t)\\
r(t) &\approx r(t + \Delta t) - \Delta t\, v(t + \Delta t) + \frac{1}{2m} \Delta t^2 F(t+\Delta t)
\end{align}

Solve for $v$:
\begin{align}
v(t+\Delta t) &\approx v(t) + \frac{1}{2m} \Delta t \big(F(t) + F(t+\Delta t)\big)
\end{align}

Complete **Velocity Verlet** integrator consists of the first and third equation.

In practice, split into three steps (calculate the velocity at the half time step):
\begin{align}
v(t+\frac{\Delta t}{2}) &= v(t) + \frac{\Delta t}{2} \frac{F(t)}{m} \\
r(t + \Delta r) &= r(t) + \Delta t\, v(t+\frac{\Delta t}{2})\\
v(t+\Delta t) &= v(t+\frac{\Delta t}{2}) + \frac{\Delta t}{2} \frac{F(t+\Delta t)}{m}
\end{align}

When writing production-level code, remember to re-use $F(t+\Delta t)$ als the "new" starting $F(t)$ in the next iteration (and don't recompute).

### Integration of planetary motion 
Gravitational potential energy:
$$
U(r) = -\frac{GMm}{r}
$$
with $r$ the distance between the two masses $m$ and $M$.

#### Central forces
$$
U(\mathbf{r}) = f(r) = f(\sqrt{\mathbf{r}\cdot\mathbf{r}})\\
\mathbf{F} = -\nabla U(\mathbf{r}) = -\frac{\partial f(r)}{\partial r} \, \frac{\mathbf{r}}{r} 
$$

#### Force of gravity
\begin{align}
\mathbf{F} &= -\frac{G m M}{r^2} \hat{\mathbf{r}}\\
\hat{\mathbf{r}} &= \frac{1}{\sqrt{x^2 + y^2}} \left(\begin{array}{c} x \\ y \end{array}\right)
\end{align}

#### Integrate simple planetary orbits 
Set $M = 1$ (one solar mass) and $m = 3.003467×10^{-6}$ (one Earth mass in solar masses) and try initial conditions

\begin{alignat}{1}
x(0) &= 1,\quad  y(0) &= 0\\
v_x(0) &= 0,\quad  v_y(0) &= 6.179
\end{alignat}

Note that we use the following units:
* length in astronomical units (1 AU = 149,597,870,700 m )
* mass in solar masses (1 $M_☉ = 1.988435×10^{30}$ kg)
* time in years (1 year = 365.25 days, 1 day = 86400 seconds)

In these units, the gravitational constant is $G = 4\pi^2$ (in SI units $G = 6.674×10^{-11}\, \text{N}\cdot\text{m}^2\cdot\text{kg}^{-2}$).


```python
M_earth = 3.003467e-6
M_sun = 1.0
```


```python
G_grav = 4*np.pi**2

def F_gravity(r, m=M_earth, M=M_sun):
    rr = np.sum(r*r)
    rhat = r/np.sqrt(rr)
    return -G_grav*m*M/rr * rhat

def U_gravity(r, m=M_earth, M=M_sun):
    return -G_grav*m*M/np.sqrt(np.sum(r*r))
```

Let's now integrate the equations of motions under gravity with the **Velocity Verlet** algorithm:


```python
def planet_orbit(r0=np.array([1, 0]), v0=np.array([0, 6.179]), mass=M_earth, dt=0.001, t_max=1):
    """2D planetary motion with velocity verlet"""
    dim = len(r0)
    assert len(v0) == dim

    nsteps = int(t_max/dt)

    r = np.zeros((nsteps, dim))
    v = np.zeros_like(r)

    r[0, :] = r0
    v[0, :] = v0

    # start force evaluation for first step
    Ft = F_gravity(r[0], m=mass)
    for i in range(nsteps-1):
        vhalf = v[i] + 0.5*dt * Ft/mass
        r[i+1, :] = r[i] + dt * vhalf
        Ftdt = F_gravity(r[i+1], m=mass)
        v[i+1] = vhalf + 0.5*dt * Ftdt/mass
        # new force becomes old force
        Ft = Ftdt
    
    return r, v
```


```python
r, v = planet_orbit(dt=0.1, t_max=10)
```


```python
rx, ry = r.T
ax = plt.subplot(1,1,1)
ax.set_aspect(1)
ax.plot(rx, ry)
```

These are not closed orbits (as we would expect from a $1/r$ potential). But it gets much besser when stepsize is reduced to 0.01 (just rerun the code with `dt = 0.01` and replot):


```python
r, v = planet_orbit(dt=0.01, t_max=10)
```


```python
rx, ry = r.T
ax = plt.subplot(1,1,1)
ax.set_aspect(1)
ax.plot(rx, ry)
```

## Velocity Verlet vs RK4: Energy conservation
Assess the stability of `rk4` and `Velocity Verlet` by checking energy conservation over longer simulation times.

The file `integrators2.py` contains almost all code that you will need.

### Implement gravity force in `integrators2.py`
Add `F_gravity` to the `integrators2.py` module. Use the new function `unitvector()`.

### Planetary orbits with `integrators2.py` 


```python
r0 = np.array([1, 0])
v0 = np.array([0, 6.179])
```


```python
from importlib import reload
reload(integrators2)
```




    <module 'integrators2_solution' from '/Volumes/Data/oliver/Documents/Teaching/ASU/CompPhys_PHY494/2018/PHY494-resources/10_ODEs/integrators2_solution.py'>



Use the new function `integrators2.integrate_newton_2d()` to integrate 2d coordinates.

#### RK4


```python
trk4, yrk4 = integrators2.integrate_newton_2d(x0=r0, v0=v0, t_max=30, mass=M_earth,
                                       h=0.01,
                                       force=integrators2.F_gravity, 
                                       integrator=integrators2.rk4)
```


```python
rxrk4, ryrk4 = yrk4[:, 0, 0], yrk4[:, 0, 1]
ax = plt.subplot(1,1,1)
ax.set_aspect(1)
ax.plot(rxrk4, ryrk4)
```


```python
integrators2.analyze_energies(trk4, yrk4, integrators2.U_gravity, m=M_earth)
```


```python
print("Energy conservation RK4 for {} steps: {}".format(
        len(trk4),
        integrators2.energy_conservation(trk4, yrk4, integrators2.U_gravity, m=M_earth)))
```

    Energy conservation RK4 for 3000 steps: 3.5366039779498594e-06


#### Euler


```python
te, ye = integrators2.integrate_newton_2d(x0=r0, v0=v0, t_max=30, mass=M_earth,
                                         h=0.01,
                            force=F_gravity, 
                            integrator=integrators2.euler)
rex, rey = ye[:, 0].T
```


```python
ax = plt.subplot(1,1,1)
ax.plot(rx, ry, label="RK4")
ax.plot(rex, rey, label="Euler")
ax.legend(loc="best")
ax.set_aspect(1)
```


```python
integrators2.analyze_energies(te, ye, integrators2.U_gravity, m=M_earth)
```


```python
print("Energy conservation Euler for {} steps: {}".format(
        len(te),
        integrators2.energy_conservation(te, ye, integrators2.U_gravity, m=M_earth)))
```

    Energy conservation Euler for 3000 steps: 0.6763156457080265


*Euler* is just awful... but we knew that already.

#### Velocity Verlet


```python
tv, yv = integrators2.integrate_newton_2d(x0=r0, v0=v0, t_max=30, mass=M_earth,
                                       h=0.01,
                                       force=F_gravity, 
                                       integrator=integrators2.velocity_verlet)
```


```python
rxv, ryv = yv[:, 0].T
ax = plt.subplot(1,1,1)
ax.set_aspect(1)
ax.plot(rxv, ryv, label="velocity Verlet")
ax.plot(rxrk4, ryrk4, label="RK4")
ax.legend(loc="best")
```


```python
integrators2.analyze_energies(tv, yv, integrators2.U_gravity, m=M_earth)
```


```python
print("Energy conservation Velocity Verlet for {} steps: {}".format(
        len(tv),
        integrators2.energy_conservation(tv, yv, integrators2.U_gravity, m=M_earth)))
```

    Energy conservation Velocity Verlet for 3000 steps: 6.37822713158452e-05


*Velocity Verlet* only has moderate accuracy, especially when compared to *RK4*.

However, let's look at energy conservation over longer times:

#### Longer time scale stability
Run RK4 and Velocity Verlet for longer.


```python
tv2, yv2 = integrators2.integrate_newton_2d(x0=r0, v0=v0, t_max=1000, mass=M_earth,
                                       h=0.01,
                                       force=F_gravity, 
                                       integrator=integrators2.velocity_verlet)
```


```python
print("Energy conservation Velocity Verlet for {} steps: {}".format(
        len(tv2),
        integrators2.energy_conservation(tv2, yv2, integrators2.U_gravity, m=M_earth)))
```

    Energy conservation Velocity Verlet for 100000 steps: 6.379584728912368e-05



```python
t4, y4 = integrators2.integrate_newton_2d(x0=r0, v0=v0, t_max=1000, mass=M_earth,
                                       h=0.01,
                                       force=F_gravity, 
                                       integrator=integrators2.rk4)
```


```python
print("Energy conservation RK4 for {} steps: {}".format(
        len(t4),
        integrators2.energy_conservation(t4, y4, integrators2.U_gravity, m=M_earth)))
```

    Energy conservation RK4 for 100000 steps: 0.00011860288209883538


Velocity Verlet shows **good long-term stability** but relative low precision. On the other hand, RK4 has high precision but the accuracy decreases over time.

* Use a *Verlet* integrator when energy conservation is important and long term stability is required (e.g. molecular dynamics simulations). It is generally recommended to use an integrator that conserves some of the inherent symmetries and structures of the governing physical equations (e.g. for Hamilton's equations of motion, time reversal symmetry and the symplectic and area-preserving structure).
* Use *RK4* for high short-term accuracy (but may be difficult to know what "short term" should mean) or when solving general differential equations.



```python

```

41a8178d4fab6a4073fa66fc996e59a7961a92e4
ipynb
Jupyter Notebook
10_ODEs/10-ODE-integrators-verlet.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2018
635a6678569406e11865c8a583a56f4a3cf2bdc4
CC-BY-4.0
10_ODEs/10-ODE-integrators-verlet.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2018
635a6678569406e11865c8a583a56f4a3cf2bdc4
CC-BY-4.0
10_ODEs/10-ODE-integrators-verlet.ipynb
ASU-CompMethodsPhysics-PHY494/PHY494-resources-2018
635a6678569406e11865c8a583a56f4a3cf2bdc4
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES