Ordinary Differential Equations

ODEs often lack analytic solution. To solve an n-th order ode numerically

convert it to a system of n first order odes
solve the n first order odes numerically

Methods to solve first order ode

Euler Method
```
y_i+1 = y_i + dy/dx * h
x_i+1 = x_i + h
```
is a first order method (error goes as $O(h^1)$)

but there is also round off error of size $\frac{\epsilon}{h}$

Minimum error $\sim\epsilon^{1/2}$ at $h\approx \epsilon^{1/2}$

Which means single precision $\implies$ $4\times 10^{-4}$

double precision $\implies$ $1\times 10^{-8}$

So most physical problems require double precision

ignores quadratic term $\implies$ poor quality
Midpoint method (2nd Order Runge-Kutta)
```
k1 = h * y'(x[i],y[i])
k2 = h * y'(x[i]+h/2,y[i]+k1/2)
y[i+1] = y[i] + k2 
```
gets error $O(h^3)$ $\implies$ second order method

You can tell from these portraits these people did not have computers. Basically they exploited women who could not get jobs as mathematicians and using them to do numerical calculations.

On Runge and Kutta, David Buscher, PhD.
Classic Runge-Kutta: 4th order

basically iterate 4 times to improve the estimate of Mean Value Theorem to an error $O(h^5)$

Comparisons

higher order methods achieve better accuracy in much fewer steps, and have lower minimum errors that low order methods never reach.
low order methods could be intrinsically unstable for certain methods, which could be tested by checking against conservation properties

Adaptive stepping

step size varied for different correlation scales along the problem

Symplectic integrators

conserve surrogate (nearby) Hamiltonian
could be particularly stable for N-body problems
not many libraries exist

Aside: Code quality

# TOY CODE TIME !!!
def method():
    """
    Docstring for a method helps users enquire documentation from command line using command doc
    what it means
    what it returns
    """
    for i in range(10):
        print('{}'.format(i)) # Good code should be self-explanatory, only subtleties should have comments
    return [0], [1], [2]

import matplotlib.pyplot as plt
fig1, ax1 = plt.subplots() # object-oriented, if use ax1.plot() instead of use plt.plot()
t, v, y = method()
ax1.plot(t,y,label='posi')
ax1.plot(t,v,label='velo')
ax1.set_ylim(down,up)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.legend()