<table>
 <tr align=left><td>
 <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>
</table>

Note:  The presentation below largely follows part II in "Finite Difference Methods for Ordinary and Partial Differential Equations" by LeVeque (SIAM, 2007).


```python
from __future__ import print_function

%matplotlib inline
import numpy
import matplotlib.pyplot as plt
```

# Numerical Solution to ODE Initial Value Problems - Part 1

Many physical, biological, and societal systems can be written as a system of ordinary differential equations (ODEs).  In the case where the initial state (value) is know the problems can be written as

$$
    \frac{\text{d}\mathbf{u}}{\text{d}t} = \mathbf{f}(t, \mathbf{u}) \quad \mathbf{u}(0) = \mathbf{u}_0
$$

where
 - $\mathbf{\!u}(t)$ is the state vector
 - $\mathbf{\!f}(t, \mathbf{\!u})$ is a vector-valued function that controls the growth of $\mathbf{u}$ with time
 - $\mathbf{\!u}(0)$ is the initial condition at time $t = 0$

#### Examples:  Simple radioactive decay
$$
    \mathbf{\!u} = [c]
$$
   
$$
    \frac{\text{d} c}{\text{d}t} = \lambda c \quad c(0) = c_0
$$
   

which has solutions of the form $c(t) = c_0 e^{\lambda t}$


```python
t = numpy.linspace(0.0, 3200., 100)
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, 1.0 * numpy.exp(decay_constant * t))
axes.plot(1600.,0.5,'ro')
axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years", fontsize=18)
axes.set_xlabel('t (years)', fontsize=16)
axes.set_ylabel('$c$', fontsize=16)
axes.set_xlim((0.0, t[-1]))
axes.set_ylim((0.0, 1.1))
axes.grid()
plt.show()
```

#### Examples:  Complex radioactive decay (or chemical system).

Chain of decays from one species to another.

$$\begin{aligned}
    \frac{\text{d} c_1}{\text{d}t} &= -\lambda_1 c_1 \\
    \frac{\text{d} c_2}{\text{d}t} &= \lambda_1 c_1 - \lambda_2 c_2 \\
    \frac{\text{d} c_3}{\text{d}t} &= \lambda_2 c_3 - \lambda_3 c_3 
\end{aligned}$$

$$\frac{\text{d} \mathbf{u}}{\text{d}t} = \frac{\text{d}}{\text{d}t}\begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix} = 
\begin{bmatrix} 
    -\lambda_1 & 0 & 0 \\
    \lambda_1 & -\lambda_2 & 0 \\
    0 & \lambda_2 & -\lambda_3
\end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix}$$

$$\frac{\text{d} \mathbf{u}}{\text{d}t} = A \mathbf{u}$$

For systems of equations like this the general solution to the ODE is the matrix exponential:

$$\mathbf{u}(t) = \mathbf{u}_0 e^{A t}$$

#### Examples:  Particle tracking in a fluid

$$\frac{\text{d} \mathbf{X}}{\text{d}t} = \mathbf{V}(t, \mathbf{X})$$

In fact all ODE IVP systems can be thought of as tracking particles through a flow field (dynamical system).  In 1-dimension the flow "manifold" we are on is fixed by the initial condition.


```python
x = numpy.linspace(0., 1., 11)
y = numpy.linspace(0., 1., 11)
x_fine = numpy.linspace(0., 1.)
y_fine = numpy.linspace(0., 1.)

X, Y = numpy.meshgrid(x,y)
X_fine, Y_fine = numpy.meshgrid(x_fine, y_fine)

pi = numpy.pi
psi = numpy.sin(pi*X_fine)*numpy.sin(pi*Y_fine)
U = pi*numpy.sin(pi*X)*numpy.cos(pi*Y)
V = -pi*numpy.cos(pi*X)*numpy.sin(pi*Y)

x0 = 0.75
y0 = 0.75
psi0 = numpy.sin(pi*x0)*numpy.sin(pi*y0)

fig = plt.figure(figsize=(6,6))
axes = fig.add_subplot(1, 1, 1)
axes.quiver(X,Y, U, V)
axes.plot(.75, 0.75,'ro')
axes.contour(X_fine, Y_fine, psi, [ psi0 ])
axes.grid()
axes.set_title("Particle tracking", fontsize=18)
axes.set_xlabel('y', fontsize=16)
axes.set_ylabel('x', fontsize=16)

plt.show()
```

#### Examples: Van der Pol Oscillator

$$y'' - \mu (1 - y^2) y' + y = 0 \quad \quad \text{with} \quad \quad  y(0) = y_0, \quad y'(0) = v_0$$
 

$$\mathbf{u} = \begin{bmatrix} y \\ y' \end{bmatrix} = \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}$$
   
$$\frac{\text{d}}{\text{d}t} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \begin{bmatrix} u_2 \\ \mu (1 - u_1^2) u_2 - u_1 \end{bmatrix} = \mathbf{f}(t, \mathbf{u})$$


```python
from scipy.integrate import solve_ivp

def f_vanderpol(t, u, mu=5):
    return numpy.array([u[1], mu * (1.0 - u[0]**2) * u[1] - u[0]])

# N = 100
N = 500
t_span = (0., 200.)
u0 = [ 1., 0. ]
f = lambda t, u : f_vanderpol(t, u, mu=50)
sol = solve_ivp(f, t_span, u0, method='BDF')
```


```python
fig = plt.figure(figsize=(16,6))
axes = fig.add_subplot(1, 2, 1)

axes.plot(sol.t, sol.y[0])
axes.set_title("Solution to Van der Pol Oscillator", fontsize=18)
axes.set_xlabel("t", fontsize=16)
axes.set_ylabel("y(t)", fontsize=16)
axes.grid()

axes = fig.add_subplot(1, 2, 2)

axes.plot(sol.y[0],sol.y[1],'r')
axes.set_title("Phase Diagram for Van der Pol Oscillator", fontsize=18)
axes.set_xlabel("y(t)", fontsize=16)
axes.set_ylabel("y'(t)", fontsize=16)
axes.grid()
plt.show()
```

## Basic Stepping Schemes

Introducing some notation to simplify things
$$\begin{aligned}
    t_0 &= 0 \\
    t_1 &= t_0 + \Delta t \\
    t_n &= t_{n-1} + \Delta t = n \Delta t + t_0 \\
    u_0 &= u(t_0) \approx U_0 \\
    u_1 &= u(t_1) \approx U_1 \\
    u_n &= u(t_n) \approx U_2 \\
\end{aligned}$$
where lower-case letters are "exact".  

Looking back at our work on numerical differentiation why not approximate the derivative as a finite difference:

$$
    \frac{u(t + \Delta t) - u(t)}{\Delta t} = f(t, u)
$$

We still need to decide how to evaluate the $f(t, u)$ term however.  

#### Relationship to quadrature

Actually it is often more instructive to work with the integral form of an ODE

$$
    \int^{t + \Delta t}_t \frac{\text{d} u}{\text{d}\tilde{\!t}} d\tilde{\!t} = \int^{t + \Delta t}_t f(t, u) d\tilde{\!t}
$$

which is equivalent to the differential form.  However using the fundamental theorem of calculus tells us that the LHS is $u(t + \Delta t) - u(t)$  so we can write the ODE as

$$ 
    u(t + \Delta t) = u(t) + \int^{t + \Delta t}_t f(\tilde{\!t}, u(\tilde{\!t})) d\tilde{\!t}
$$




#### The Geometric picture

$$
    \frac{\text{d} u}{\text{d} t} = f(t,u)
$$
says $f$ is a slope field that is tangent to any solution $u(t)$ that passes through $u_0$ at $t=t_0$


```python
t = numpy.linspace(0., 4800, 11)
y = numpy.linspace(0., 1.2, 11)
T, Y = numpy.meshgrid(t,y)
dt = numpy.ones(T.shape)
dy = -dt*Y

tp = numpy.linspace(0., 4800, 100)
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.quiver(T,Y, dt,dy)
axes.plot(tp,numpy.exp(decay_constant*tp))
axes.plot(0.,1.,'ro')
axes.grid()
axes.set_title("Direction Set, $u' = - \lambda u$, $u(0)=1$", fontsize=18)
axes.set_xlabel('t (years)', fontsize=16)
axes.set_ylabel('u', fontsize=16)

axes.set_ylim((-.1,1.2))
plt.show()


```

Alternatively 

$$ 
    u(t + \Delta t) = u(t) + \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t})) d\tilde{\!t}
$$
Implies, that $u$ at some time $\Delta t$ in the future,  is $u(t)$ plus a *number* 
$$
    K = \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t}) )d\tilde{\!t}
$$
which is a definite *line integral* (along an unknown solution).  As we will see many of the algorithms for integrating ODE's can be related to simple quadrature rules for approximating $K$

## Single-Step Multi-Stage Schemes

The integral form of an ODE initial value problem can be written

$$ 
    u(t + \Delta t) = u(t) + \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t})) d\tilde{\!t}
$$

Which says that our solution $u$, if it exists  at some time $\Delta t$ in the future,  is $u(t)$ plus a *number* 
$$
    K = \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t}) )d\tilde{\!t}
$$
which is a definite *line integral* (along an unknown solution). 

An important class of ODE solvers are called *Single Step, Multi-stage schemes* which can be most easily understood as extensions of the  Newton-Cotes quadrature schemes for approximating $K$ (plus an error term that will scale as $\Delta t^p$)

### The Geometric picture



```python
t = numpy.linspace(0., 4800, 11)
y = numpy.linspace(0., 1.2, 11)
T, Y = numpy.meshgrid(t,y)
dt = numpy.ones(T.shape)
dy = -dt*Y

tK = 2000.
uK = numpy.exp(decay_constant*tK)
K = uK -1.
tp = numpy.linspace(0., 4800, 100)
tk = numpy.linspace(0., tK, 100)

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.quiver(T,Y, dt,dy)
axes.plot(tp,numpy.exp(decay_constant*tp))
axes.plot(0.,1.,'ro')
axes.plot(tk,numpy.exp(decay_constant*tk),'r--')
axes.plot(tK, uK, 'ro')
axes.plot([0.,0.], [1., uK], 'r--')
axes.text(10., 0.72, '$K$', fontsize=24, color='red')
axes.plot([0.,tK],[uK, uK], 'r--')
axes.text(900., uK - .1, '$\Delta t$', fontsize=24, color='red')
axes.text(-10, 1., '$U_0$', fontsize=24, color='blue')
axes.text(tK+10, uK, '$U_1$', fontsize=24, color='blue')


axes.grid()
axes.set_title("Direction Set, $u' = - \lambda u$, $u(0)=1$", fontsize=18)
axes.set_xlabel('t (years)', fontsize=16)
axes.set_ylabel('u', fontsize=16)

axes.set_ylim((-.1,1.2))
plt.show()


```

#### Forward Euler scheme
For example, if we approximate $K$ with a left-sided quadrature rule
$$
    K = \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t})) d\tilde{\!t} \approx \Delta t f(t, u(t))
$$

then our first ODE algorithm can be written
$$
u(t + \Delta t) =  u(t) + \Delta t f(t, u(t))
$$


or in terms of our discrete approximation $U$
$$
\begin{align}
    K_1 &= \Delta t f(t_n, U_n)\\
    U_{n+1} &= U_n + K_1\\
\end{align}
$$

which is known as the *forward Euler method*.  In essence we are approximating the derivative with the value of the function at the point we are at $t_n$.


```python
t = numpy.linspace(0.0, 1.6e3, 100)
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0

# Euler step
dt = 1e3
u_np = c_0 + dt * (decay_constant * c_0)
```


```python
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, c_0 * numpy.exp(decay_constant * t), label="True Solution")
axes.plot(0., 1., 'ro')
axes.text(0., 1.01, '$U_0$', fontsize=16)
axes.plot((0.0, dt), (c_0, u_np), 'k')
axes.plot((dt, dt), (u_np, c_0 * numpy.exp(decay_constant * dt)), 'k--')
axes.plot((0.0, 0.0), (c_0, u_np), 'k--')
axes.text(10., 0.75, '$K_1$', fontsize=16)
axes.plot((0.0, dt), (u_np, u_np), 'k--')
axes.text(400, u_np - 0.05, '$\Delta t$', fontsize=16)

axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.05))
axes.grid()
plt.show()
```


```python
# Implement Forward Euler
def euler(f, t_span, u0, N):
    """ simple implementation of constant step-size forward euler method
        This doc string should have so much more in it
    """
    t = numpy.linspace(t_span[0], t_span[1],N)
    u = numpy.empty(t.shape)
    u[0] = u0
    delta_t = t[1] - t[0]
    for (n, t_n) in enumerate(t[:-1]):
        K1 = delta_t * f(t_n, u[n])
        u[n + 1] = u[n] + K1        
    return t, u
```


```python
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

t_span = [0.0, 1.6e3]
u0 = 1.
N = 20
t_euler, u_euler = euler(f, t_span, u0, N)

```


```python
t_exact = numpy.linspace(0.0, 1.6e3, 100)
u_exact = c_0 * numpy.exp(decay_constant * t_exact)

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t_euler, u_euler, 'or', label="Euler")
axes.plot(t_exact, u_exact, 'k--', label="True Solution")

axes.set_title("Forward Euler")
axes.set_xlabel("t (years)")
axes.set_xlabel("$c(t)$")
axes.set_ylim((0.4,1.1))
axes.grid()
axes.legend()
plt.show()
```

### Backward's Euler

Similar to forward Euler is the *backward Euler* method which,   uses a right-rectangle rule to estimate $K$ given $f$ at a future time. i.e. 

$$
    K\approx \Delta t f(t_{n+1}, U_{n+1})
$$

However, the update scheme now becomes
$$
    U_{n+1} = U_n + \Delta t f(t_{n+1}, U_{n+1}).
$$

which requires a (usually non-linear) solve for $U_{n+1}$. Schemes where the function $f$ is evaluated at the unknown time are called *implicit methods*.

For some cases we can solve the equation by hand.  For instance in the case of our example problem we have:

$$
    U_{n+1} = U_n + \Delta t f(t_{n+1}, U_{n+1}) = U_n + \Delta t (\lambda U_{n+1})
$$

which can be solved for $U_{n+1}$ to find

$$\begin{aligned}
    U_{n+1} &= U_n + \Delta t (\lambda U_{n+1}) \\
    U_{n+1} \left[ 1 - \Delta t \lambda \right ] &= U_n \\
    U_{n+1} &= \frac{U_n}{1 - \Delta t \lambda}
\end{aligned}$$


```python
t = numpy.linspace(0.0, 1.6e3, 100)
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, c_0 * numpy.exp(decay_constant * t), label="True Solution")

# Plot Backwards Euler step
dt = 1e3
u_np = c_0 + dt * (decay_constant * c_0 * numpy.exp(decay_constant * dt))
axes.plot((0.0, dt), (c_0, u_np), 'k')
axes.plot(dt, u_np, 'ro')
axes.text(dt+ 10., u_np, '$U_1$', fontsize=16)
axes.plot((0.0, 0.0), (c_0, u_np), 'k--')
axes.plot((0.0, dt), (u_np, u_np), 'k--')
axes.text(400, u_np - 0.05, '$\Delta t$', fontsize=16)
axes.text(10., 0.85, '$K_1$', fontsize=16)

axes.grid()
axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.05))
plt.show()
```


```python
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

t_exact = numpy.linspace(0.0, 1.6e3, 100)
u_exact = c_0 * numpy.exp(decay_constant * t_exact)

# Implement backwards Euler
t_backwards = numpy.linspace(0.0, 1.6e3, 10)
delta_t = t_backwards[1] - t_backwards[0]
u_backwards = numpy.empty(t_backwards.shape)
u_backwards[0] = c_0
for n in range(0, t_backwards.shape[0] - 1):
    u_backwards[n + 1] = u_backwards[n] / (1.0 - decay_constant * delta_t)

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t_backwards, u_backwards, 'or', label="Backwards Euler")
axes.plot(t_exact, u_exact, 'k--', label="True Solution")
axes.grid()
axes.set_title("Backwards Euler")
axes.set_xlabel("t (years)")
axes.set_xlabel("$c(t)$")
axes.set_ylim((0.4,1.1))
axes.legend()
plt.show()
```

It's also useful to be able to do this in the case of systems of ODEs.  Let $f(U) = A U$, then

$$\begin{aligned}
    U_{n+1} &= U_n + \Delta t (A U_{n+1}) \\
    U_{n+1} \left [ I - \Delta t A \right ] &= U_n \\
    U_{n+1} &= \left [ I - \Delta t A \right]^{-1} U_n
\end{aligned}$$

In general however we are often not able to do this with arbitrary $f$.

Another simple implicit method is based on quadrature using the trapezoidal method.  The scheme is
$$
    \frac{U_{n+1} - U_{n}}{\Delta t} = \frac{1}{2} (f(U_n) + f(U_{n+1}))
$$

In this case what is the update scheme?

$$\begin{aligned}
    U_{n+1} &= U_{n} + \frac{\Delta t}{2} (f(U_n) + f(U_{n+1})) \\
    U_{n+1} &= U_{n} + \frac{\Delta t}{2} (\lambda U_n + \lambda U_{n+1}) \\
    U_{n+1} \left[1 - \frac{\Delta t \lambda}{2}  \right] &= U_{n} \left[1 + \frac{\Delta t \lambda}{2} \right] \\
    U_{n+1} &= U_{n} \frac{1 + \frac{\Delta t \lambda}{2}}{1 - \frac{\Delta t \lambda}{2}} \\
\end{aligned}$$


```python
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0
t_exact = numpy.linspace(0.0, 1.6e3, 100)
u_exact = c_0 * numpy.exp(decay_constant * t_exact)

# Implement trapezoidal method
t = numpy.linspace(0.0, 1.6e3, 10)
delta_t = t[1] - t[0]
u = numpy.empty(t.shape)
u[0] = c_0
integration_constant = (1.0 + decay_constant * delta_t / 2.0) / (1.0 - decay_constant * delta_t / 2.0)
for n in range(t.shape[0] - 1):
    u[n + 1] = u[n] * integration_constant

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, u, 'or', label="Trapezoidal")
axes.plot(t_exact, u_exact, 'k--', label="True Solution")
axes.grid()

axes.set_title("Trapezoidal")
axes.set_xlabel("t (years)")
axes.set_xlabel("$c(t)$")
axes.set_ylim((0.4,1.1))
axes.legend()
plt.show()
```

## Error Analysis of ODE Methods

At this point it is also helpful to introduce more notation to distinguish between the true solution to the ODE $u(t_n)$ and the approximated value which we will denote $U_n$.

**Definition:** We define the *truncation error* of a scheme by replacing the $U_n$ with the true solution $u(t_n)$ in the finite difference formula and looking at the difference from the exact solution.

For example we will use the difference form of forward Euler
$$
    \frac{U_{n+1} - U_n}{\Delta t} = f(t_n)
$$
and define the truncation error as
$$
    T(t, u; \Delta t) = \frac{u(t_{n+1}) - u(t_n)}{\Delta t} - f(t_n, u(t_n)).
$$

**Definition:** A method is called *consistent* if 
$$
    \lim_{\Delta t \rightarrow 0} T(t, u; \Delta t) = 0.
$$

**Definition:** We say that a method is *order* $p$ accurate if

$$
    \lVert T(t, u; \Delta t) \rVert \leq C \Delta t^p
$$

uniformally on $t \in [0, T]$.  This can also be written as $T(t, u; \Delta t) = \mathcal{O}(\Delta t^p)$.  Note that a method is consistent if $p > 0$.

### Error Analysis of Forward Euler

We can analyze the error and convergence order of forward Euler by considering the Taylor series centered at $t_n$:

$$
    u(t) = u(t_n) + (t - t_n) u'(t_n) + \frac{u''(t_n)}{2} (t - t_n)^2 + \mathcal{O}((t-t_n)^3)
$$

Evaluating this series at $t_{n+1}$ gives

$$\begin{aligned}
    u(t_{n+1}) &= u(t_n) + (t_{n+1} - t_n) u'(t_n) + \frac{u''(t_n)}{2} (t_{n+1} - t_n)^2 + \mathcal{O}((t_{n+1}-t_n)^3)\\
    &=u_n + \Delta t f(t_n, u_n) + \frac{u''(t_n)}{2} \Delta t^2 + \mathcal{O}(\Delta t^3)
\end{aligned}$$

From the definition of truncation error we can use our Taylor series expression and find the truncation error.  Take the finite difference form of forward Euler

$$
    \frac{U_{n+1} - U_n}{\Delta t} = f(t_n)
$$

and replacing the derivative formulation with $u(t_n)$ to find

$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{u(t_{n+1}) - u(t_n)}{\Delta t} - f(t_n) \\
    &= \frac{u(t_{n+1}) - u(t_n)}{\Delta t} - u'(t_n).
\end{aligned}$$

From here we use the Taylor series centered at $t_n$ and evaluated at $t_{n+1}$ to find

$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{u(t_{n+1}) - u(t_n)}{\Delta t} - u'(t_n) \\
    &= \frac{1}{\Delta t} \left[ u(t_n) + u'(t_n) (t - t_n) + \frac{u''(t_n)}{2} (t - t_n)^2 + \mathcal{O}((t-t_n)^3) - u(t_n) \right] - u'(t_n) \\
    &=  u'(t_n) + \frac{u''(t_n)}{2} \Delta t + \mathcal{O}(\Delta t^2) - u'(t_n) \\
    &= \frac{u''(t_n)}{2} \Delta t + \mathcal{O}(\Delta t^2).
\end{aligned}$$

This implies that forward Euler is first order accurate and therefore consistent.

Another equivalent definition of the truncation error uses the form
$$
    U_{n+1} = u(t_n) + \Delta t f(t_n)
$$
and the definition
$$
    T(t, u; \Delta t) = \frac{1}{\Delta t} \left [ U_{n+1} - u(t_{n+1}) \right]
$$

to find
$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{1}{\Delta t} [U_{n+1} - u(t + \Delta t)] \\
    &= \frac{1}{\Delta t} \left[ \underbrace{u_n + \Delta t f(t_n, u_n)}_{U_{n+1}} - \underbrace{\left( u_n + \Delta t f(t_n, u_n) + \frac{u''(t_n)}{2} \Delta t^2 + \mathcal{O}(\Delta t^3) \right )}_{u(t_{n+1})}\right ] \\
    &= \frac{1}{\Delta t} \left[ - \frac{u''(t_n)}{2} \Delta t^2 - \mathcal{O}(\Delta t^3) \right ] \\
    &= - \frac{u''(t_n)}{2} \Delta t - \mathcal{O}(\Delta t^2)
\end{aligned}$$

#### Truncation Error vs Step Error

Sometimes we will also consider the "Step Error"  which is the error that is introduced over one step

$$
    E_h = | U_{n+1} - u_{n+1} |
$$

which in general will be of order $O(\Delta t^{p+1})$ if the Truncation error is $O(\Delta t^p)$.  So for Forward (or Backward's) Euler the step error 

$$
    E_h = O(\Delta t^2)
$$

The step error can be very useful in *adaptive stepping* schemes

## Runge-Kutta Methods

One way to derive higher-order ODE solvers is to use higher order quadrature schemes that sample the function at a number of  intermediate stages to provide a more accurate estimate of $K$.  These are not *multi-step* methods as they still only require information from the current time step but they raise the order of accuracy by adding *stages*.  These types of methods are called **Runge-Kutta** methods.

### Example:  Two-stage Runge-Kutta Methods

The basic idea behind the simplest of the Runge-Kutta methods is to approximate $K$ using a mid-point scheme (which should be 2nd order accurate.  Unforunately, we don't know the value of the mid-point.  However we can use an Euler step of size $\Delta_t/2$ to estimate the mid-point.  

We can write the algorithm as 

$$\begin{aligned}
    K_1 &= \Delta t f(U_n, t_n) \\
    K_2 &= \Delta t f(U_n + K_1/2, t_n + \Delta/2, )\\
    U_{n+1} &= U_n + K_2 \\    
\end{aligned}$$

Where we now evaluate the function in two stages $K_1$ and $K_2$.

or for an autonomous ODE
$$
    U_{n+1} = U_n + \Delta t f(U_n + \frac{1}{2} \Delta t f(U_n))
$$


```python
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

# RK2 step
dt = 1e3
U0 = 1.0
K1 = dt * f(0., U0)
Y1 = U0 + K1/2
K2 = dt * f(dt/2., Y1)
U1 = U0 + K2


fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, U0 * numpy.exp(decay_constant * t), label="True Solution")
axes.plot(0., U0, 'ro')
axes.text(0., U0+.01, '$U_0$', fontsize=16)
axes.plot((0.0, dt), (U0, U0 + K1), 'k--')
axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k')

axes.plot(dt/2., U0 + K1/2, 'ro')
axes.plot((0.0, 0.0), (U0, Y1), 'k--')
axes.text(10., 0.85, '$\\frac{K_1}{2}$', fontsize=18)
axes.plot((0.0, dt/2), (Y1, Y1), 'k--')
axes.text(250, Y1 - 0.05, '$\\frac{\Delta t}{2}$', fontsize=18)

axes.plot(dt, U1, 'go')
axes.plot((0.0, 0.0), (U0, Y1), 'k--')
axes.text(10., 0.85, '$\\frac{K_1}{2}$', fontsize=18)
axes.plot((0.0, dt/2), (Y1, Y1), 'k--')
axes.text(250, Y1 - 0.05, '$\\frac{\Delta t}{2}$', fontsize=18)

axes.plot(dt, U1, 'go')
axes.plot((0., dt), (U0, U1), 'k')
axes.text(dt+20, U1, '$U_1$', fontsize=18)
#axes.plot((0.0, 0.0), (U0, U1), 'g--')
#axes.plot((0.0, dt), (U1, U1), 'g--')



axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.05))
axes.grid()
plt.show()
```

#### Error analysis RK2

The truncation error can be computed similarly to how we did so before but we do need to figure out how to compute the derivative inside of the function.  Note that due to 
$$
    f(u(t_n)) = u'(t_n)
$$ 
that differentiating this leads to 
$$
    f'(u(t_n)) u'(t_n) = u''(t_n)
$$ 
leading to
$$\begin{aligned}
    f\left(u(t_n) + \frac{1}{2} \Delta t f(u(t_n)) \right ) &= f\left(u(t_n) +\frac{1}{2} \Delta t u'(t_n) \right ) \\
    &= f(u(t_n)) + \frac{1}{2} \Delta t u'(t_n) f'(u(t_n)) + \frac{1}{8} \Delta t^2 (u'(t_n))^2 f''(u(t_n)) + \mathcal{O}(\Delta t^3) \\
    &=u'(t_n) + \frac{1}{2} \Delta t u''(t_n) + \mathcal{O}(\Delta t^2)
\end{aligned}$$

Going back to the truncation error we have
$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{1}{\Delta t} \left[u_n + \Delta t f\left(u_n + \frac{1}{2} \Delta t f(u_n)\right) - \left(u_n + \Delta t f(t_n, u_n) + \frac{u''(t_n)}{2} \Delta t^2 + \mathcal{O}(\Delta t^3) \right ) \right] \\
    &=\frac{1}{\Delta t} \left[\Delta t u'(t_n) + \frac{1}{2} \Delta t^2 u''(t_n) + \mathcal{O}(\Delta t^3) - \Delta t u'(t_n) - \frac{u''(t_n)}{2} \Delta t^2 + \mathcal{O}(\Delta t^3) \right] \\
    &= \mathcal{O}(\Delta t^2)
\end{aligned}$$

so this method is second order accurate.

### Example:  4-stage Runge-Kutta Method

If RK2 is related to a Mid-point quadrature scheme,  then the classic 4-stage, 4th order Runge-Kutta scheme should be reminiscent of Simpson's Quadrature rule.  It requires 4 samples of $f(t,u)$ at the beginning of the step, two-samples in the middle and one at the end, then a linear combination of those samples

$$\begin{aligned}
    K_1 &= \Delta t f(t_n, U_n) \\
    K_2 &= \Delta t f(t_n + \Delta t/2, U_n + K_1/2) \\
    K_3 &= \Delta t f(t_n + \Delta t/2, U_n + K_2/2) \\
    K_4 &= \Delta t f(t_n + \Delta t, U_n + K_3) \\
        & \\
    U_{n+1} &= U_n + \frac{1}{6} \left [K_1 + 2(K_2 + K_3)  + K_4) \right ] 
\end{aligned}$$

With truncation error $T = O(\Delta t^4)$


```python
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

# RK4 step
dt = 1e3
U0 = 1.0
K1 = dt * f(0., U0)
K2 = dt * f(dt/2., U0 + K1/2)
K3 = dt * f(dt/2., U0 + K2/2)
K4 = dt * f(dt, U0 + K3)

U1 = U0 + 1./6. *( K1 + 2 * (K2 + K3) + K4)


fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, U0 * numpy.exp(decay_constant * t), label="True Solution")
axes.plot(0., U0, 'ro')
axes.text(0.-20, U0-.04, '$K_1$', color='red',fontsize=16)

axes.text(0., U0+.01, '$U_0$', fontsize=16)
axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k--')
axes.plot(dt/2., U0 + K1/2, 'ro')
axes.text(dt/2-20, U0 + K1/2-.04, '$K_2$', color='red',fontsize=16)


axes.plot((0.0, dt/2.), (U0, U0 + K2/2.), 'k--')
axes.plot(dt/2., U0 + K2/2, 'ro')
axes.text(dt/2-20, U0 + K2/2+.02, '$K_3$', color='red',fontsize=16)

axes.plot((0.0, dt), (U0, U0 + K3), 'k--')
axes.plot(dt, U0 + K3, 'ro')
axes.text(dt-20, U0 + K3-.04, '$K_4$', color='red',fontsize=16)

axes.plot(dt, U1, 'go')
#axes.plot((0., dt), (U0, U1), 'k')
axes.text(dt+20, U1, '$U_1$', fontsize=18)



axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.05))
axes.grid()
plt.show()
```


```python
def RK2(f, t_span, u0, N):
    """ implement constant step size 2 stage Runge-Kutta Method RK2"""
    
    t = numpy.linspace(t_span[0], t_span[1], N)
    delta_t = t[1] - t[0]
    u = numpy.empty(t.shape)
    u[0] = u0 
    for (n, t_n) in enumerate(t[1:]):
        K_1 = delta_t * f(t_n, u[n])
        K_2 = delta_t * f(t_n + delta_t/2., u[n] + K_1/2.)
        u[n+1] = u[n] + K_2
    return t, u

def improved_euler(f, t_span, u0, N):
    """ implement constant step size 2 stage Improved Euler Method trapezoidal rule"""
    
    t = numpy.linspace(t_span[0], t_span[1], N)
    delta_t = t[1] - t[0]
    u = numpy.empty(t.shape)
    u[0] = u0
    for (n, t_n) in enumerate(t[1:]):
        K_1 = delta_t * f(t_n, u[n])
        K_2 = delta_t * f(t_n + delta_t, u[n] + K_1)
        u[n+1] = u[n] + 0.5 * (K_1 + K_2)
    return t, u 
```


```python
def RK4(f, t_span, u0, N):
    """ implement constant step size 4 stage Runge-Kutta Method RK4"""
    
    t = numpy.linspace(t_span[0], t_span[1], N)
    delta_t = t[1] - t[0]
    u = numpy.empty(t.shape)
    u[0] = u0
    for (n, t_n) in enumerate(t[1:]):
        K_1 = delta_t * f(t_n, u[n])
        K_2 = delta_t * f(t_n + delta_t/2., u[n] + K_1/2.)
        K_3 = delta_t * f(t_n + delta_t/2., u[n] + K_2/2.)
        K_4 = delta_t * f(t_n + delta_t, u[n] + K_3)
        u[n+1] = u[n] + 1./6. * (K_1 + 2.*( K_2 + K_3) + K_4)
    return t, u
```


```python
# Implement and compare the two-stage and 4-stage Runge-Kutta methods
f = lambda t, u: -u

N = 20
t_span = [ 0., 5.0 ]
u0 = 1.

t_exact = numpy.linspace(t_span[0], t_span[1], 100)
u_exact = u0*numpy.exp(-t_exact)
t_euler, u_euler = euler(f, t_span, u0, N)
t_ieuler, u_ieuler = improved_euler(f, t_span, u0, N)
t_RK2, u_RK2 = RK2(f, t_span, u0, N)
t_RK4, u_RK4 = RK4(f, t_span, u0, N)
```


```python
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t_exact,u_exact,'k',label='exact')
axes.plot(t_euler, u_euler, 'ro', label='euler')
axes.plot(t_ieuler, u_ieuler, 'co', label='improved euler')
axes.plot(t_RK2, u_RK2, 'go', label='RK2')
axes.plot(t_RK4, u_RK4, 'bo', label='RK4')

axes.grid()
axes.set_xlabel('t', fontsize=16)
axes.set_ylabel('u', fontsize=16)
axes.legend(loc='best')
plt.show()

```

### Convergence of Single Step Multi-Stage schemes

All of the above schemes are consistent and have truncation errors $T\propto\Delta t^p$


```python
N = numpy.array([ 2**n for n in range(4,10)])
err_euler = numpy.zeros(len(N))
err_ieuler = numpy.zeros(len(N))
err_RK2 = numpy.zeros(len(N))
err_RK4 = numpy.zeros(len(N))

t_span = [ 0., 4.]
dt = t_span[1]/N

u0 = 1. 
u_exact = u0*numpy.exp(-t_span[1])

for i, n in enumerate(N):
    t, u_euler = euler(f, t_span, u0, n)
    err_euler[i] = numpy.abs(u_euler[-1] - u_exact)
    t, u_ieuler = improved_euler(f, t_span, u0, n)
    err_ieuler[i] = numpy.abs(u_ieuler[-1] - u_exact)
    t, u_RK2 = RK2(f, t_span, u0, n)
    err_RK2[i] = numpy.abs(u_RK2[-1] - u_exact)
    t, u_RK4 = RK4(f, t_span, u0, n)
    err_RK4[i] = numpy.abs(u_RK4[-1] - u_exact)
    
err_fit = lambda dt, p: numpy.exp(p[1])*dt**p[0]
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

# Euler
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_euler[2:]),1)
line = axes.loglog(dt, err_euler, 'o', label='euler, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

# Improved Euler
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_ieuler[2:]),1)
line = axes.loglog(dt, err_ieuler, 'o', label='improved euler, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

# RK2
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK2[2:]),1)
line = axes.loglog(dt, err_RK2, 'o', label='rk2, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

#RK4
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK4[2:]),1)
line = axes.loglog(dt, err_RK4, 'o', label='rk4, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())


axes.grid()
axes.set_xlabel('$\Delta t$', fontsize=16)
axes.set_ylabel('$Error$', fontsize=16)
axes.set_title('Convergence: Single Step Schemes', fontsize=18)
axes.legend(loc='best', fontsize=14)
plt.show()
```

## Summary: single-Step Multi-Stage Schemes

The integral form of an ODE initial value problem can be written

$$ 
    u(t + \Delta t) = u(t) + \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t})) d\tilde{\!t}
$$

Which says that our solution $u$, if it exists  at some time $\Delta t$ in the future,  is $u(t)$ plus a *number* 
$$
    K = \int^{t + \Delta t}_t f(\tilde{t}, u(\tilde{t}) )d\tilde{\!t}
$$
which is a definite *line integral* (along an unknown solution). 

#### Single Step, Multi-stage schemes

are most easily understood as extensions of the  Newton-Cotes quadrature schemes for approximating $K$ (plus an error term that will scale as $\Delta t^p$)

**Explicit Schemes**
<table width="80%">
    <tr align="center"><th>Name</th> <th align="center">Stages</th> <th align="center">"Quadrature"</th><th align="center">$$T$$</th></tr>
     <tr align="center"><td>Euler</td> <td align="center">1</td> <td align="center">Left-Rectangle</td><td align="center">$$O(\Delta t)$$</td></tr>
    <tr align="center"><td>Improved Euler</td> <td align="center">2</td> <td align="center">Trapezoidal</td><td align="center">$$O(\Delta t^2)$$</td></tr>
    <tr align="center"><td>RK2</td> <td align="center">2</td> <td align="center">Mid-Point</td><td align="center">$$O(\Delta t^2)$$</td></tr>
    <tr align="center"><td>RK4</td> <td align="center">4</td> <td align="center">Simpson</td><td align="center">$$O(\Delta t^4)$$</td></tr>
</table>

**Implicit Schemes**
<table width="80%">
    <tr align="center"><th>Name</th> <th align="center">Stages</th> <th align="center">"Quadrature"</th><th align="center">$$T$$</th></tr>
     <tr align="center"><td>Backwards-Euler</td> <td align="center">1</td> <td align="center">Right-Rectangle</td><td align="center">$$O(\Delta t)$$</td></tr>
    <tr align="center"><td>Trapezoidal</td> <td align="center">2</td> <td align="center">Trapezoidal</td><td align="center">$$O(\Delta t^2)$$</td></tr>
</table>


```python
N = numpy.array([ 2**n for n in range(4,10)])
err_euler = numpy.zeros(len(N))
err_ieuler = numpy.zeros(len(N))
err_RK2 = numpy.zeros(len(N))
err_RK4 = numpy.zeros(len(N))

t_span = [ 0., 4.]
dt = t_span[1]/N

u0 = 1. 
u_exact = u0*numpy.exp(-t_span[1])

for i, n in enumerate(N):
    t, u_euler = euler(f, t_span, u0, n)
    err_euler[i] = numpy.abs(u_euler[-1] - u_exact)
    t, u_ieuler = improved_euler(f, t_span, u0, n)
    err_ieuler[i] = numpy.abs(u_ieuler[-1] - u_exact)
    t, u_RK2 = RK2(f, t_span, u0, n)
    err_RK2[i] = numpy.abs(u_RK2[-1] - u_exact)
    t, u_RK4 = RK4(f, t_span, u0, n)
    err_RK4[i] = numpy.abs(u_RK4[-1] - u_exact)
    
err_fit = lambda dt, p: numpy.exp(p[1])*dt**p[0]
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

# Euler
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_euler[2:]),1)
line = axes.loglog(dt, err_euler, 'o', label='euler, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

# Improved Euler
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_ieuler[2:]),1)
line = axes.loglog(dt, err_ieuler, 'o', label='improved euler, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

# RK2
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK2[2:]),1)
line = axes.loglog(dt, err_RK2, 'o', label='rk2, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())

#RK4
p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK4[2:]),1)
line = axes.loglog(dt, err_RK4, 'o', label='rk4, p={:3.2f}'.format(p[0]))
axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())


axes.grid()
axes.set_xlabel('$\Delta t$', fontsize=16)
axes.set_ylabel('$Error$', fontsize=16)
axes.set_title('Convergence: Single Step Schemes', fontsize=18)
axes.legend(loc='best', fontsize=14)
plt.show()
```

## Adaptive Time Stepping

#### Why should we care about all of these schemes and their errors?

* Even though we know the formal error.  It is with respect to a true solution we don't know.
* In itself, the error estimates don't tell us how to choose a time step $\Delta t$ to keep the solution within a given tolerance
* However,  in combination, we can use multiple methods to control the error and provide **Adaptive** time stepping

#### Example:  Compare 1 step of Euler to one step of RK2



```python
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

# RK2 step
dt = 1e3
U0 = 1.0
K1 = dt * f(0., U0)
Y1 = U0 + K1/2
K2 = dt * f(dt/2., Y1)
U1 = U0 + K2

t = numpy.linspace(U0, 1600.)
u_true = U0 * numpy.exp(decay_constant * t)

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, u_true, label="True Solution")
axes.plot(0., U0, 'ro')
axes.text(0., U0+.01, '$U_0$', fontsize=16)
axes.plot((0.0, dt), (U0, U0 + K1), 'k')
#axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k')

#Euler step
axes.plot(dt, U0 + K1, 'ro')
#axes.plot((0.0, 0.0), (U0, Y1), 'k--')
axes.text(dt + 10., U0 + K1, '$U_{euler}$', fontsize=18)

axes.plot(dt, U1, 'go')
#axes.plot((0.0, 0.0), (U0, Y1), 'k--')
#axes.text(10., 0.85, '$\\frac{K_1}{2}$', fontsize=18)
#axes.plot((0.0, dt/2), (Y1, Y1), 'k--')
#axes.text(250, Y1 - 0.05, '$\\frac{\Delta t}{2}$', fontsize=18)

# RK2 Step
axes.plot(dt, U1, 'go')
axes.plot((0., dt), (U0, U1), 'k')
axes.text(dt+20, U1, '$U_{RK2}$', fontsize=18)
#axes.plot((0.0, 0.0), (U0, U1), 'g--')
#axes.plot((0.0, dt), (U1, U1), 'g--')

# difference
axes.plot((dt, dt), (U1, U0+K1),'k--')
axes.text(dt+40, 0.5*(U1 + U0+K1), '$\Delta\propto\Delta t^2$', fontsize=18)


axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.05))
axes.legend(loc='best')
axes.grid()
plt.show()
```

#### Relative Truncation Error 

If we consider the *Step Error* for each of our schemes, we know that

$$
\begin{align}
    u_{n+1} &= U^{euler}_{n+1} + O(\Delta t^2)\\
    u_{n+1} &= U^{RK2}_{n+1} + O(\Delta t^3)\\
\end{align}
$$

Therefore we can compute the *relative truncation error* as

$$
    \Delta = | U^{euler}_{n+1} - U^{RK2}_{n+1} | = O(\Delta t^{?})
$$

* $\Delta$ is Computable!
* has a known dependence on step-size $\Delta t$

### Adaptive Time Stepping

Given the relative truncation error and its scaling with $\Delta t$, we can now use this to choose a single good time step.

#### Example:

Suppose we wanted our relative truncation error to be small relative to the solution or zero, we could set 

$$
    \Delta_{target} = \mathtt{rtol}\,U^{RK2}_{n+1} + \mathtt{atol}
$$

where $\mathtt{rtol}$ and $\mathtt{atol}$ are relative and absolute tolerances (and we assume that $U^{RK2}_{n+1}$ is a reasonably good estimate of the true solution)


But our measured relative truncation error, $\Delta$ depends on  the step size we just took i.e

$$
    \Delta_{measured} \propto \Delta t_n^2
$$

### Adaptive Time Stepping

If we take the ratio of both relationships we get

$$
    \frac{\Delta_{target}}{\Delta_{measured}} = \left[\frac{\Delta t_{target}}{\Delta t_{n}}\right]^2
$$

or rearranging, our target step size is

$$
    \Delta t_{target} = \Delta t_{n}\left[\frac{\Delta_{target}}{\Delta_{measured}}\right]^{\frac{1}{2}}
$$

which tells us how to grow or shrink our time step to maintain accuracy.


In general,  if we have two methods with different step errors such that 

$$
    \Delta \propto \Delta t^p
$$

then our adaptive stepper will look like

$$
    \Delta t_{target} = \Delta t_{n}\left[\frac{\Delta_{target}}{\Delta_{measured}}\right]^{1/p}
$$

This leads to all sorts of adaptive schemes most are included in standard libraries.  

### Embedded Runge-Kutta Schemes

There are in fact a whole family of **Embedded RK** schemes which are $N$ stage schemes but can combine the $N$ function evaluations in two different ways to produce methods with different error estimates.  

A popular one is **RK45** (available in `SciPy`) which is based on the [Dormand-Prince 5(4)](https://doi.org/10.1016/0771-050X(80)90013-3) pair which uses 6 function evaluations per step to produce a 4th order and 5th order scheme.  The 4th order scheme controls the time step, and the 5th order scheme actually is the solution.


```python
from scipy.integrate import solve_ivp

def f_vanderpol(t, u, mu=5):
    return numpy.array([u[1], mu * (1.0 - u[0]**2) * u[1] - u[0]])

t_span = (0., 50.)
u0 = [ 1., 0. ]
f = lambda t, u : f_vanderpol(t, u, mu=20) 
sol = solve_ivp(f, t_span, u0, method='RK45',rtol=1.e-3,atol=1.e-6)
```


```python
fig = plt.figure(figsize=(16,6))
axes = fig.add_subplot(1, 2, 1)

axes.plot(sol.t, sol.y[0],'o-')
axes.set_title("Solution to Van der Pol Oscillator", fontsize=18)
axes.set_xlabel("t", fontsize=16)
axes.set_ylabel("y(t)", fontsize=16)
axes.grid()

axes = fig.add_subplot(1, 2, 2)
delta_t = sol.t[1:] - sol.t[:-1]
axes.plot(sol.t[:-1], delta_t)
axes.grid()
axes.set_xlabel('$t$', fontsize=16)
axes.set_ylabel('$\Delta t$', fontsize=16)
axes.set_title('Timesteps, N = {}'.format(len(sol.t)), fontsize=18)

plt.show()
```

## Taylor Series Methods

A **Taylor series method** can be derived by direct substitution of the right-hand-side function $f(t, u)$ and its appropriate derivatives into the Taylor series expansion for $u(t_{n+1})$.  For a $p$th order method we would look at the Taylor series up to that order and replace all the derivatives of $u$ with derivatives of $f$ instead.  

For the general case we have
$$\begin{align*}
    u(t_{n+1}) = u(t_n) + \Delta t u'(t_n) + \frac{\Delta t^2}{2} u''(t_n) + \frac{\Delta t^3}{6} u'''(t_n) + \cdots + \frac{\Delta t^p}{p!} u^{(p)}(t_n)
\end{align*}$$
which contains derivatives of $u$ up to $p$th order.  

We then replace these derivatives with the appropriate derivative of $f$ which will always be one less than the derivative of $u$ (due to the original ODE)

$$
    u^{(p)}(t_n) = f^{(p-1)}(t_n, u(t_n))
$$

leading to the method

$$
\begin{align}
    u(t_{n+1}) &= u(t_n) + \Delta t f(t_n, u(t_n)) + \frac{\Delta t^2}{2} f'(t_n, u(t_n)) \\
    &+ \frac{\Delta t^3}{6} f''(t_n, u(t_n)) + \cdots + \frac{\Delta t^p}{p!} f^{(p-1)}(t_n, u(t_n)).
\end{align}
$$

### 2nd Order Taylor Series Method

We want terms up to second order so we need to take the derivative of $u' = f(t, u)$ once to find $u'' = f'(t, u)$ and therefore
$$\begin{align*}
    u(t_{n+1}) &= u(t_n) + \Delta t u'(t_n) + \frac{\Delta t^2}{2} u''(t_n) \\
    &=u(t_n) + \Delta t f(t_n, u(t_n)) + \frac{\Delta t^2}{2} f'(t_n, u(t_n)) ~~~ \text{or} \\
    U_{n+1} &= U_n + \Delta t f(t_n, U_n) + \frac{\Delta t^2}{2} f'(t_n, U_n).
\end{align*}$$

With Step error $O(\Delta t^3)$ and truncation error $T$, $O(\Delta t^2)$

### Example

Let's use our simplest problem $u'(t) = \lambda u$ with $f=\lambda u$.  Therefore

$$\begin{align*}
    f(t,u) &= \lambda u\\
    f'(t,u) &= \lambda u' = \lambda f = \lambda^2 u\\
    f''(t,u) &= \lambda^2 u' = \lambda^2 f = \lambda^3 u
\end{align*}$$

so a third order scheme would look like
$$
\begin{align}
    U(t_{n+1}) &= U(t_n)\left[ 1 + \lambda\Delta t  + \frac{(\lambda\Delta t)^2}{2} + \frac{(\lambda\Delta t)^3}{6}\right]+ O(\Delta t^4)
\end{align}
$$


```python
def Taylor_3_flambda_u(lamda, t_span, u0, N):
    """ implement constant step size  3rd order Taylor Series method for f(t,u) = \lambda u"""
    
    t = numpy.linspace(t_span[0], t_span[1], N)
    lambda_dt = lamda*(t[1] - t[0])
    u = numpy.empty(t.shape)
    u[0] = u0 
    for (n, t_n) in enumerate(t[1:]):
        u[n+1] = u[n] * ( 1. + lambda_dt  + (lambda_dt**2)/2.  + (lambda_dt**3)/6.)
    return t, u

```


```python
lam = -1.
t_span = [0., 5.]
u0 = 1.

f = lambda t,u : -u
t_exact = numpy.linspace(t_span[0], t_span[1], 100)
u_exact = u0*numpy.exp(-t_exact)

N = 20.
t_taylor, u_taylor = Taylor_3_flambda_u(lam, t_span, u0, N)
t_euler, u_euler = euler(f, t_span, u0, N)
```


```python
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t_exact,u_exact,'k',label='exact')
axes.plot(t_euler, u_euler, 'ro', label='euler')
axes.plot(t_taylor, u_taylor, 'bo', label='Taylor3')

axes.grid()
axes.set_xlabel('t', fontsize=16)
axes.set_ylabel('u', fontsize=16)
axes.legend(loc='best')
plt.show()

```

### Some Drawbacks

**Taylor Series methods**
 - require differentiating the given equation which can be cumbersome and difficult to implement
 - require a new routine for every $f$
 
**General one-step/multi-stage methods**
  - higher order methods often require a large number of evaluations of $f$ per time step

## Linear Multi-Step Methods

**Multi-step methods**  are ODE methods that 
 - require only *one* new function evaluation per time step  to work.  
 - reuse values and function evaluations at some number of previous time steps
 
**Disadvantages over single step methods**

 - Methods are not self-starting, i.e. they require other methods to find the initial values
 - Difficult to adapt. The time step  Δ𝑡  in one-step methods can be changed at any time while multi-step methods this is much more complex

### Simplest example:  The leap-frog method

The leap-frog method is similar to Euler's method in that it uses the information from two-previous time steps to advance the problem.  We can write the problem as a centered first-derivative centered   at time $t_{n+1}, U_{n+1}$ i.e. 

$$\frac{U_{n+2} - U_{n}}{2\Delta t} = f(t_{n+1}, U_{n+1})$$

or

$$
    U_{n+2} = U_{n} + 2\Delta t\, f(t_{n+1}, U_{n+1})
$$ 

this method is known as the leap-frog method


```python
t = numpy.linspace(0.0, 1.6e3, 100)
c_0 = 1.0
decay_constant = -numpy.log(2.0) / 1600.0

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t, c_0 * numpy.exp(decay_constant * t), label="True Solution")

# Plot Leap-Frog step
dt = 1e3
u1 = c_0 * numpy.exp(decay_constant * dt / 2.0)
u_np = c_0 + dt * (decay_constant * u1)
axes.plot((0.0, dt), (c_0, u_np), 'k')
axes.plot((dt, dt), (u_np, c_0 * numpy.exp(decay_constant * dt)), 'k--')
axes.plot((0.0, 0.0), (c_0, u_np), 'k--')
axes.plot((0.0, dt), (u_np, u_np), 'k--')
axes.text(400, u_np - 0.05, '$\Delta t$', fontsize=16)
axes.plot([0., dt/2, dt], [ c_0, u1, u_np],'ro')
axes.set_title("Radioactive Decay with $t_{1/2} = 1600$ years")
axes.set_xlabel('t (years)')
axes.set_ylabel('$c$')
axes.grid()
axes.set_xlim(-1e2, 1.6e3)
axes.set_ylim((0.5,1.0))
plt.show()
```


```python
def leap_frog(f, t_span, u0, N, start=RK2):
    """ calculate fixed step with leap-frog iterator with a single step starter
    """
    t = numpy.linspace(t_span[0], t_span[1], N)
    delta_t = t[1] - t[0]
    u = numpy.zeros(t.shape)
    u[0] = u0
    # use a single-step multi-stage method to start
    
    t_start, u_start = start(f, (t[0],t[1]), u0, 2)
    u[1] = u_start[-1]
    for (n, t_np) in enumerate(t[1:-1]):
        u[n+2] = u[n] +  2 *delta_t * f(t_np, u[n+1]) 
    return t, u
```


```python
u0 = 1.0
t_span = (0., 1600.) 
N = 7

# Stable example
decay_constant = -numpy.log(2.0) / 1600.0
f = lambda t, u: decay_constant * u

t_exact = numpy.linspace(t_span[0], t_span[1], N)
u_exact = u0 * numpy.exp( decay_constant * t_exact)

t_leapfrog, u_leapfrog = leap_frog(f, t_span, u0, N, start=RK2)
```


```python
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)
axes.plot(t_leapfrog, u_leapfrog, 'or-', label="Leap-Frog")
axes.plot(t_exact, u_exact, 'k--', label="True Solution")
axes.grid()
axes.set_title("Leap-Frog", fontsize=18)
axes.set_xlabel("t (years)", fontsize=16)
axes.set_xlabel("$c(t)$", fontsize=16)
axes.legend(loc='best', fontsize=14)
plt.show()
```

### Error Analysis of Leap-Frog Method

To easily analyze this method we will expand the Taylor series around time $t_{n+1}$ to yield
$$\begin{aligned}
    u(t_{n+2}) &= u_{n+1} + \Delta t f(t_{n+1},u_{n+1}) + \Delta t^2 \frac{u''(t_{n+1})}{2}  + \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4)
\end{aligned}$$

We need one more expansion however due to leap-frog.  Recall that leap-frog has the form
$$
    U_{n+2} = U_{n} + 2 \Delta t f(t_{n+1}, U_{n+1}).
$$
To handle the $U_{n}$ term we need to write this with relation to $u(t_{n+1})$.  Again we use the Taylor series
$$
    u(t_n) = u_{n+1} - \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  - \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4)
$$

$$\begin{aligned}
    u(t_{n+2}) &= u_{n+1} + \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  + \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4)\\
    u(t_{n}) &= u_{n+1} - \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  - \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4)
\end{aligned}$$

Plugging these into our definition of the truncation error along with the leap-frog method definition leads to

$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{1}{\Delta t} \left [\underbrace{U_{n} + 2 \Delta t f_{n+1}}_{U_{n+2}} - \underbrace{\left(u_{n+1} + \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  + \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4) \right )}_{u(t_{n+2})} \right ] \\
    &=\frac{1}{\Delta t} \left [ \underbrace{ \left(u_{n+1} - \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  - \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4)\right)}_{U_{n}} + 2\Delta t f_n - \underbrace{\left(u_{n+1} + \Delta t f_{n+1} + \Delta t^2 \frac{u''(t_{n+1})}{2}  + \Delta t^3 \frac{u'''(t_{n+1})}{6} + \mathcal{O}(\Delta t^4) \right )}_{u(t_{n+2})} \right ] \\
    &=\frac{1}{\Delta t} \left [- \Delta t^3 \frac{u'''(t_n)}{3} + \mathcal{O}(\Delta t^4) \right ] \\ 
        &=- \Delta t^2 \frac{u'''(t_n)}{3} + \mathcal{O}(\Delta t^3)
\end{aligned}$$

Therefore the method is second order accurate and is consistent theoretically.  In practice it's a bit more complicated than that.


```python
# Compare accuracy between Euler, RK2 and Leap-Frog
f = lambda t, u: -u
u_exact = lambda t: numpy.exp(-t)
u_0 = 1.0

t_span = (0.0, 10.0)
num_steps = [2**n for n in range(4,11)]
```


```python
delta_t = numpy.empty(len(num_steps))
error_euler = numpy.empty(len(num_steps))
error_RK2 = numpy.empty(len(num_steps))
error_leapfrog = numpy.empty(len(num_steps))

for (i, N) in enumerate(num_steps):
    t = numpy.linspace(t_span[0], t_span[1], N)
    tt, u_euler = euler(f, t_span, u_0, N )
    tt, u_rk2 = RK2(f, t_span, u_0, N)
    tt, u_leapfrog = leap_frog(f, t_span, u_0, N, start=euler)
    
    delta_t[i] = t[1] - t[0]
        
    # Compute error for each
    error_euler[i] = numpy.linalg.norm(delta_t[i] * (u_euler - u_exact(t)), ord=1)
    error_RK2[i] = numpy.linalg.norm(delta_t[i] * (u_rk2 - u_exact(t)), ord=1)
    error_leapfrog[i] = numpy.linalg.norm(delta_t[i] * (u_leapfrog - u_exact(t)), ord=1)
    
# Plot error vs. delta_t
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))
axes.loglog(delta_t, error_euler, 'bo', label='Forward Euler')
axes.loglog(delta_t, error_RK2, 'ro', label='RK2')
axes.loglog(delta_t, error_leapfrog, 'go', label="Leap-Frog")

axes.loglog(delta_t, order_C(delta_t[2], error_euler[2], 1.0) * delta_t**1.0, '--b')
axes.loglog(delta_t, order_C(delta_t[2], error_RK2[2], 2.0) * delta_t**2.0, '--r')
axes.loglog(delta_t, order_C(delta_t[2], error_leapfrog[2], 2.0) * delta_t**2.0, '--r')

axes.grid()
axes.legend(loc=2, fontsize=14)
axes.set_title("Comparison of Errors", fontsize=18)
axes.set_xlabel("$\Delta t$",fontsize=16)
axes.set_ylabel("$|U(t_f) - u(t_f)|$", fontsize=16)

plt.show()
```

### Look at the errors for Leap-Frog

They're actually quite large...If you make a quick plot of u_leapfrog vs $t$ you'll see what's happening (and is a good example of an issue we will need to address in future lectures)


```python
N= 1000
t_leapfrog, u_leapfrog = leap_frog(f, t_span, u_0, N, start=euler)
## Your plotting code here
plt.figure()
plt.plot(t_leapfrog, u_leapfrog)
plt.grid()
plt.show()

```

### General Linear Multi-Step Methods

Leap-frog is perhaps the simplest of multi-step methods but all linear multi-step methods can be written as the linear combination of past, present and future solutions:

Want to get from time $U_n$ to time $U_{n+r}$ by taking steps $j$:

$$
    \sum^r_{j=0} \alpha_j U_{n+j} = \Delta t \sum^r_{j=0} \beta_j f(U_{n+j}, t_{n+j})
$$
If $\beta_r = 0$ then the method is explicit (only requires previous time steps).  Note that the coefficients are not unique as we can multiply both sides by a constant.  In practice a normalization of $\alpha_r = 1$ is used.

For example:  our Leap-frog method can be written using $r=2$, $\alpha = \begin{bmatrix} -1 & 0 & 1\\
\end{bmatrix}$, $\beta = \begin{bmatrix} 0 & 2 & 0 \\ \end{bmatrix}$

#### Example: Adams Methods

$$
    U_{n+r} = U_{n+r-1} + \Delta t \sum^r_{j=0} \beta_j f(U_{n+j}).
$$
All these methods have $\alpha_r = 1$, $\alpha_{r-1} = -1$ and $\alpha_j=0$ for $j < r - 1$.

### Adams-Bashforth Methods
The **Adams-Bashforth** methods are explicit solvers that maximize the order of accuracy given a number of steps $r$.  This is accomplished by looking at the Taylor series and picking the coefficients $\beta_j$ to eliminate as many terms in the Taylor series as possible. AB is initialized with forward Euler.
$$\begin{aligned}
    \text{1-step:} & ~ & U_{n+1} &= U_n +\Delta t f(U_n) \\
    \text{2-step:} & ~ & U_{n+2} &= U_{n+1} + \frac{\Delta t}{2} (-f(U_n) + 3 f(U_{n+1})) \\
    \text{3-step:} & ~ & U_{n+3} &= U_{n+2} + \frac{\Delta t}{12} (5 f(U_n) - 16 f(U_{n+1}) + 23 f(U_{n+2})) \\
    \text{4-step:} & ~ & U_{n+4} &= U_{n+3} + \frac{\Delta t}{24} (-9 f(U_n) + 37 f(U_{n+1}) -59 f(U_{n+2}) + 55 f(U_{n+3}))
\end{aligned}$$


```python
def AB2(f, t_span, u0, N, start=RK2):
    """ calculate fixed step Adams-Bashforth 2-step method with a single step starter
        reuses previous function evaluations
    """
    t = numpy.linspace(t_span[0], t_span[1], N)
    delta_t = t[1] - t[0]
    u = numpy.zeros(t.shape)
    
    u[0] = u0
    # use a single-step multi-stage method to start
    t_start, u_start = start(f, (t[0],t[1]), u0, 2)
    u[1] = u_start[-1]
    
    # set initial function evaluations
    fn = f(t[0], u[0])
    fnp = f(t[1], u[1])
    for (n, t_np) in enumerate(t[2:]):
        u[n+2] = u[n + 1] + delta_t / 2.0 * (-fn + 3.0 * fnp)
        fn = fnp
        fnp = f(t_np, u[n+2])
    return t, u
```


```python
# Use 2-step Adams-Bashforth to compute solution
f = lambda t, u: -u
t_span = (0., 10.)
u0 = 1.0

N = 20
t, u_ab2 = AB2(f, t_span, u0, N, start=RK2)
```


```python
t_exact = numpy.linspace(t_span[0], t_span[1], 100)
u_exact = numpy.exp(-t_exact)

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

axes.plot(t_exact, u_exact, 'k', label="True")
axes.plot(t, u_ab2, 'ro', label="2-step A-B")

axes.set_title("Adams-Bashforth Method", fontsize=18)
axes.set_xlabel("t", fontsize=16)
axes.set_ylabel("u(t)",fontsize=16)
axes.legend(loc=1, fontsize=14)
axes.grid()

plt.show()
```

### Adams-Moulton Methods
The **Adams-Moulton** methods are the implicit versions of the Adams-Bashforth methods.  Since this gives one additional parameter to use $\beta_r$ these methods are generally one order of accuracy greater than their counterparts. AM is initialized with trapezoidal scheme.
$$\begin{aligned}
    \text{1-step:} & ~ & U_{n+1} &= U_n + \frac{\Delta t}{2} (f(U_n) + f(U_{n+1})) \\
    \text{2-step:} & ~ & U_{n+2} &= U_{n+1} + \frac{\Delta t}{12} (-f(U_n) + 8f(U_{n+1}) + 5f(U_{n+2})) \\
    \text{3-step:} & ~ & U_{n+3} &= U_{n+2} + \frac{\Delta t}{24} (f(U_n) - 5f(U_{n+1}) + 19f(U_{n+2}) + 9f(U_{n+3})) \\
    \text{4-step:} & ~ & U_{n+4} &= U_{n+3} + \frac{\Delta t}{720}(-19 f(U_n) + 106 f(U_{n+1}) -264 f(U_{n+2}) + 646 f(U_{n+3}) + 251 f(U_{n+4}))
\end{aligned}$$


```python
# Use 2-step Adams-Moulton to compute solution
# u' = - decay u
decay_constant = 1.0
f = lambda t, u: -decay_constant * u

t_exact = numpy.linspace(0.0, 10.0, 100)
u_exact = numpy.exp(-t_exact)

N = 20
# N = 10
# N = 5
t = numpy.linspace(0, 10.0, N)
delta_t = t[1] - t[0]
U = numpy.empty(t.shape)
U[0] = 1.0
U[1] = U[0] + 0.5 * delta_t * f(t[0], U[0])
U[1] = U[0] + delta_t * f(t[0], U[1])    
integration_constant = 1.0 / (1.0 + 5.0 * decay_constant * delta_t / 12.0)
for n in range(t.shape[0] - 2):
    U[n+2] = (U[n+1] + decay_constant * delta_t / 12.0 * (U[n] - 8.0 * U[n+1])) * integration_constant

fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

axes.plot(t_exact, u_exact, 'k', label="True")
axes.plot(t, U, 'ro', label="2-step A-M")

axes.set_title("Adams-Moulton Method ($f=-u$)", fontsize=18)
axes.set_xlabel("t", fontsize=16)
axes.set_ylabel("u(t)", fontsize=16)
axes.legend(loc=1, fontsize=14)
axes.grid()
plt.show() 
```

### Truncation Error for Multi-Step Methods

We can again find the truncation error in general for linear multi-step methods:
$$\begin{aligned}
    T(t, u; \Delta t) &= \frac{1}{\Delta t} \left [\sum^r_{j=0} \alpha_j u_{n+j} - \Delta t \sum^r_{j=0} \beta_j f(u_{n+j}, t_{n+j}) \right ]
\end{aligned}$$

Using the general expansion and evaluation of the Taylor series about $t_n$ we have
$$\begin{aligned}
    u(t_{n+j}) &= u(t_n) + j \Delta t u'(t_n) + \frac{1}{2} (j \Delta t)^2 u''(t_n) + \mathcal{O}(\Delta t^3) \\
    u'(t_{n+j}) &= u'(t_n) + j \Delta t u''(t_n) + \frac{1}{2} (j \Delta t)^2 u'''(t_n) + \mathcal{O}(\Delta t^3)
\end{aligned}$$

collecting terms of order $u^{(p)}$

$$
\begin{aligned}
    T(t, u; \Delta t) &= \frac{1}{\Delta t}\left( \sum^r_{j=0} \alpha_j\right) u(t_n) + \left(\sum^r_{j=0} (j\alpha_j - \beta_j)\right) u'(t_n) + \Delta t \left(\sum^r_{j=0} \left (\frac{1}{2}j^2 \alpha_j - j \beta_j \right) \right) u''(t_n) \\
& \quad \quad + \cdots + \Delta t^{q - 1} \left (\sum^r_{j=0} \left(\frac{1}{q!} j^q \alpha_j - \frac{1}{(q-1)!} j^{q-1} \beta_j \right) \right) u^{(q)}(t_n) + \cdots
\end{aligned}$$

The method is *consistent* if the first two terms of the expansion vanish, i.e. $\sum^r_{j=0} \alpha_j = 0$ and $\sum^r_{j=0} j \alpha_j = \sum^r_{j=0} \beta_j$.


```python
# Compare accuracy between RK-2, AB-2 and AM-2, RK-4
f = lambda t, u: -u
u_exact = lambda t: numpy.exp(-t)

t_f = 10.0
t_span = (0.0, t_f)

num_steps = [2**n for n in range(4,10)]
delta_t = numpy.empty(len(num_steps))
error_rk = numpy.empty(len(num_steps))
error_rk4 = numpy.empty(len(num_steps))
error_ab = numpy.empty(len(num_steps))
error_am = numpy.empty(len(num_steps))

for (i, N) in enumerate(num_steps):
    t = numpy.linspace(0, t_f, N)
    delta_t[i] = t[1] - t[0]
        
    # Compute RK2
    tt, U_rk  = RK2(f, t_span, u0, N)
    
    # Compute RK4
    tt, U_rk4  = RK4(f, t_span, u0, N)
        
    # Compute Adams-Bashforth 2-stage
    tt, U_ab  = AB2(f, t_span, u0, N)
    
    # Compute Adama-Moulton 2-stage
    U_am = numpy.empty(t.shape)
    U_am[:2] = U_rk[:2]
    decay_constant = 1.0
    integration_constant = 1.0 / (1.0 + 5.0 * decay_constant * delta_t[i] / 12.0)
    for n in range(t.shape[0] - 2):
        U_am[n+2] = (U_am[n+1] + decay_constant * delta_t[i] / 12.0 * (U_am[n] - 8.0 * U_am[n+1])) * integration_constant
        
    # Compute error for each
    error_rk[i] = numpy.linalg.norm(delta_t[i] * (U_rk - u_exact(t)), ord=1)
    error_rk4[i] = numpy.linalg.norm(delta_t[i] * (U_rk4 - u_exact(t)), ord=1)
    error_ab[i] = numpy.linalg.norm(delta_t[i] * (U_ab - u_exact(t)), ord=1)
    error_am[i] = numpy.linalg.norm(delta_t[i] * (U_am - u_exact(t)), ord=1)
    
# Plot error vs. delta_t
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

axes.loglog(delta_t, error_rk, 'ko', label='RK-2')
axes.loglog(delta_t, error_ab, 'bo', label='AB-2')
axes.loglog(delta_t, error_am, 'go', label="AM-2")
axes.loglog(delta_t, error_rk4, 'co', label='RK-4')


order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))
axes.loglog(delta_t, order_C(delta_t[0], error_ab[0], 1.0) * delta_t**1.0, '--r')
axes.loglog(delta_t, order_C(delta_t[1], error_ab[1], 2.0) * delta_t**2.0, '--b')
axes.loglog(delta_t, order_C(delta_t[1], error_am[1], 3.0) * delta_t**3.0, '--g')
axes.loglog(delta_t, order_C(delta_t[1], error_rk4[1], 4.0) * delta_t**4.0, '--c')


axes.legend(loc=4, fontsize=14)
axes.set_title("Comparison of Errors",fontsize=18)
axes.set_xlabel("$\Delta t$",fontsize=16)
axes.set_ylabel("$|U(t) - u(t)|$", fontsize=16)
axes.grid()

plt.show()
```

### Predictor-Corrector Methods

One way to simplify the Adams-Moulton methods so that implicit evaluations are not needed is by estimating the required implicit function evaluations with an explicit method.  These are often called **predictor-corrector** methods as the explicit method provides a *prediction* of what the solution might be and the now explicit *corrector* step works to make that estimate more accurate.

#### Example: One-Step Adams-Bashforth-Moulton

Use the One-step Adams-Bashforth method to predict the value of $U_{n+1}$ and then use the Adams-Moulton method to correct that value:
$$\begin{aligned}
    \hat{U}_{n+1} &= U_n + \Delta t f(U_n) \\
    U_{n+1} &= U_n + \frac{1}{2} \Delta t \left[f(U_n) + f(\hat{U}_{n+1}) \right]
\end{aligned}$$
leading to a second order accurate method.  Note this algorithm is identical to \_\_\_\_\_\_\_\_\_\_\_\__________?


```python
# One-step Adams-Bashforth-Moulton
f = lambda t, u: -u

t_exact = numpy.linspace(0.0, 10.0, 100)
u_exact = numpy.exp(-t_exact)

N = 50
t = numpy.linspace(0, 10.0, N)
delta_t = t[1] - t[0]
U = numpy.empty(t.shape)

U[0] = 1.0
for n in range(t.shape[0] - 1):
    U[n+1] = U[n] + delta_t * f(t[n], U[n])
    U[n+1] = U[n] + 0.5 * delta_t * (f(t[n], U[n]) + f(t[n+1], U[n+1]))
    
t_ie, u_ieuler = improved_euler(f, (0.0, 10.), 1., N)
fig = plt.figure(figsize=(8,6))
axes = fig.add_subplot(1, 1, 1)

axes.plot(t_exact, u_exact, 'k', label="True")
axes.plot(t, U, 'ro', label="2-step A-B")
axes.plot(t_ie, u_ieuler, 'bx', label="Improved Euler")


axes.set_title("Adams-Bashforth-Moulton P/C Method", fontsize=18)
axes.set_xlabel("t", fontsize=18)
axes.set_ylabel("u(t)", fontsize=18)
axes.legend(loc='best', fontsize=14)
axes.grid()

plt.show()
```

a8edce472413a97747b6ad4ee1695d6e81faab93
ipynb
Jupyter Notebook
09_ODE_ivp_part1.ipynb
janice-yang/intro-numerical-methods
bed376b6fe2098caba46e61c1ec53c65dba96b71
CC-BY-4.0
09_ODE_ivp_part1.ipynb
janice-yang/intro-numerical-methods
bed376b6fe2098caba46e61c1ec53c65dba96b71
CC-BY-4.0
09_ODE_ivp_part1.ipynb
janice-yang/intro-numerical-methods
bed376b6fe2098caba46e61c1ec53c65dba96b71
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES