On Schemes for Exponential Decay

class: center, middle, inverse

# On Schemes for Exponential Decay

###**Hans Petter Langtangen** at Center for Biomedical Computing, Simula Research Laboratory and Department of Informatics, University of Oslo

### Sep 24, 2015

.center[<img src="fig/CN_logo.png" width=300>]

---
class: inverse

## Goal

The primary goal of this demo talk is to demonstrate how to write
talks with [DocOnce](https://github.com/hplgit/doconce)
and get them rendered in numerous HTML formats.

---
class: inverse

# Problem setting and methods

.center[<img src="fig/method.png" width=600>]

---
class: inverse

## We aim to solve the (almost) simplest possible differential equation problem

$$
\begin{equation}
u'(t) = -au(t)
\label{ode}
\end{equation}
$$

$$
\begin{equation}  
u(0)  = I
\label{initial:value}
\end{equation}
$$

Here,

* \$ t\in (0,T] \$
 * \$ a \$, \$ I \$, and \$ T \$ are prescribed parameters
 * \$ u(t) \$ is the unknown function
 * The ODE \eqref{ode} has the initial condition \eqref{initial:value}

.center[<img src="fig/teacher2.jpg" width=250>]
---
class: inverse

## The ODE problem is solved by a finite difference scheme

* Mesh in time: \$ 0= t_0< t_1 \cdots < t_N=T \$
 * Assume constant \$ \Delta t = t_{n}-t_{n-1} \$
 * \$ u^n \$: numerical approx to the exact solution at \$ t_n \$

The \$ \theta \$ rule,

$$
u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n,
\quad n=0,1,\ldots,N-1
$$

contains the [Forward Euler](http://en.wikipedia.org/wiki/Forward_Euler_method) (\$ \theta=0 \$),
the [Backward Euler](http://en.wikipedia.org/wiki/Backward_Euler_method) (\$ \theta=1 \$),
and the [Crank-Nicolson](http://en.wikipedia.org/wiki/Crank-Nicolson) (\$ \theta=0.5 \$)
schemes.
---
class: inverse

## The Forward Euler scheme explained

---
class: inverse

## Implementation

*Implementation in a Python function:*

```python
def solver(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T]; step: dt."""
    dt = float(dt)           # avoid integer division
    N = int(round(T/dt))     # no of time intervals
    T = N*dt                 # adjust T to fit time step dt
    u = zeros(N+1)           # array of u[n] values
    t = linspace(0, T, N+1)  # time mesh

u[0] = I                 # assign initial condition
    for n in range(0, N):    # n=0,1,...,N-1
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
    return u, t
```

---
class: inverse

## How to use the solver function

*A complete main program.*

```python
# Set problem parameters
I = 1.2
a = 0.2
T = 8
dt = 0.25
theta = 0.5

from solver import solver, exact_solution
u, t = solver(I, a, T, dt, theta)

import matplotlib.pyplot as plt
plt.plot(t, u, t, exact_solution)
plt.legend(['numerical', 'exact'])
plt.show()
```

---
class: inverse

# Results

.center[<img src="fig/results.jpg" width=600>]

---
class: inverse

## The Crank-Nicolson method shows oscillatory behavior for not sufficiently small time steps, while the solution should be monotone

.center[<img src="fig/CN.png" width=600>]

---
class: inverse

## The artifacts can be explained by some theory

Exact solution of the scheme:

$$
u^n = A^n,\quad A = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}\thinspace .
$$

Key results:

* Stability: \$ |A| < 1 \$
 * No oscillations: \$ A>0 \$
 * \$ \Delta t < 1/a \$ for Forward Euler (\$ \theta=0 \$)
 * \$ \Delta t < 2/a \$ for Crank-Nicolson (\$ \theta=1/2 \$)

*Concluding remarks:* 
Only the Backward Euler scheme is guaranteed to always give
qualitatively correct results.