$$\newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\acc}{\boldsymbol{a}} \newcommand{\rpos}{\boldsymbol{r}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\q}{\boldsymbol{q}} \newcommand{\g}{\boldsymbol{g}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\stepzero}{*} \newcommand{\stephalf}{***} \newcommand{\stepone}{**} \newcommand{\baspsi}{\psi} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integer}{\mathbb{Z}}$$

# Finite difference operator notation

\begin{align} u'(t_n) &\approx \lbrack D_tu\rbrack^n = \frac{u^{n+\half} - u^{n-\half}}{\Delta t} \tag{6.1}\\ u'(t_n) &\approx \lbrack D_{2t}u\rbrack^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t} \tag{6.2}\\ u'(t_n) &= \lbrack D_t^-u\rbrack^n = \frac{u^{n} - u^{n-1}}{\Delta t} \tag{6.3}\\ u'(t_n) &\approx \lbrack D_t^+u\rbrack^n = \frac{u^{n+1} - u^{n}}{\Delta t} \tag{6.4}\\ u'(t_{n+\theta}) &= \lbrack \bar D_tu\rbrack^{n+\theta} = \frac{u^{n+1} - u^{n}}{\Delta t} \tag{6.5}\\ u'(t_n) &\approx \lbrack D_t^{2-}u\rbrack^n = \frac{3u^{n} - 4u^{n-1} + u^{n-2}}{2\Delta t} \tag{6.6}\\ u''(t_n) &\approx \lbrack D_tD_t u\rbrack^n = \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2} \tag{6.7}\\ u(t_{n+\half}) &\approx \lbrack \overline{u}^{t}\rbrack^{n+\half} = \half(u^{n+1} + u^n) \tag{6.8}\\ u(t_{n+\half})^2 &\approx \lbrack \overline{u^2}^{t,g}\rbrack^{n+\half} = u^{n+1}u^n \tag{6.9}\\ u(t_{n+\half}) &\approx \lbrack \overline{u}^{t,h}\rbrack^{n+\half} = \frac{2}{\frac{1}{u^{n+1}} + \frac{1}{u^n}} \tag{6.10}\\ u(t_{n+\theta}) &\approx \lbrack \overline{u}^{t,\theta}\rbrack^{n+\theta} = \theta u^{n+1} + (1-\theta)u^n , \tag{6.11}\\ &\qquad t_{n+\theta}=\theta t_{n+1} + (1-\theta)t_{n-1} \tag{6.12} \end{align}

Some may wonder why $$\theta$$ is absent on the right-hand side of (6.5). The fraction is an approximation to the derivative at the point $$t_{n+\theta}=\theta t_{n+1} + (1-theta) t_{n}$$.

# Truncation errors of finite difference approximations

\begin{align} \uex'(t_n) &= [D_t\uex]^n + R^n = \frac{\uex^{n+\half} - \uex^{n-\half}}{\Delta t} +R^n\nonumber,\\ R^n &= -\frac{1}{24}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4) \tag{6.13}\\ \uex'(t_n) &= [D_{2t}\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n-1}}{2\Delta t} + R^n\nonumber,\\ R^n &= -\frac{1}{6}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4) \tag{6.14}\\ \uex'(t_n) &= [D_t^-\uex]^n +R^n = \frac{\uex^{n} - \uex^{n-1}}{\Delta t} +R^n\nonumber,\\ R^n &= -\half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2) \tag{6.15}\\ \uex'(t_n) &= [D_t^+\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n}}{\Delta t} +R^n\nonumber,\\ R^n &= \half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2) \tag{6.16}\\ \uex'(t_{n+\theta}) &= [\bar D_t\uex]^{n+\theta} +R^{n+\theta} = \frac{\uex^{n+1} - \uex^{n}}{\Delta t} +R^{n+\theta}\nonumber,\\ R^{n+\theta} &= -\half(1-2\theta)\uex''(t_{n+\theta})\Delta t + \frac{1}{6}((1 - \theta)^3 - \theta^3)\uex'''(t_{n+\theta})\Delta t^2 + \nonumber\\ &\quad {\cal O}(\Delta t^3) \tag{6.17}\\ \uex'(t_n) &= [D_t^{2-}\uex]^n +R^n = \frac{3\uex^{n} - 4\uex^{n-1} + \uex^{n-2}}{2\Delta t} +R^n\nonumber,\\ R^n &= \frac{1}{3}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^3) \tag{6.18}\\ \uex''(t_n) &= [D_tD_t \uex]^n +R^n = \frac{\uex^{n+1} - 2\uex^{n} + \uex^{n-1}}{\Delta t^2} +R^n\nonumber,\\ R^n &= -\frac{1}{12}\uex''''(t_n)\Delta t^2 + {\cal O}(\Delta t^4) \tag{6.19} \end{align} \begin{align} \uex(t_{n+\theta}) &= [\overline{\uex}^{t,\theta}]^{n+\theta} +R^{n+\theta} = \theta \uex^{n+1} + (1-\theta)\uex^n +R^{n+\theta},\nonumber\\ R^{n+\theta} &= -\half\uex''(t_{n+\theta})\Delta t^2\theta (1-\theta) + {\cal O}(\Delta t^3) \tp \tag{6.20} \end{align}

# Finite differences of exponential functions

### Complex exponentials

Let $$u^n = \exp{(i\omega n\Delta t)} = e^{i\omega t_n}$$. \begin{align} [D_tD_t u]^n &= u^n \frac{2}{\Delta t}(\cos \omega\Delta t - 1) = -\frac{4}{\Delta t}\sin^2\left(\frac{\omega\Delta t}{2}\right), \tag{6.21}\\ [D_t^+ u]^n &= u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1), \tag{6.22}\\ [D_t^- u]^n &= u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}), \tag{6.23}\\ [D_t u]^n &= u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)}, \tag{6.24}\\ [D_{2t} u]^n &= u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)} \tag{6.25} \tp \end{align}

### Real exponentials

Let $$u^n = \exp{(\omega n\Delta t)} = e^{\omega t_n}$$. \begin{align} [D_tD_t u]^n &= u^n \frac{2}{\Delta t}(\cos \omega\Delta t - 1) = -\frac{4}{\Delta t}\sin^2\left(\frac{\omega\Delta t}{2}\right), \tag{6.26}\\ [D_t^+ u]^n &= u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1), \tag{6.27}\\ [D_t^- u]^n &= u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}), \tag{6.28}\\ [D_t u]^n &= u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)}, \tag{6.29}\\ [D_{2t} u]^n &= u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)} \tag{6.30} \tp \end{align}

# Finite differences of $$t^n$$

The following results are useful when checking if a polynomial term in a solution fulfills the discrete equation for the numerical method. \begin{align} \lbrack D_t^+ t\rbrack^n = 1, \tag{6.31}\\ \lbrack D_t^- t\rbrack^n = 1, \tag{6.32}\\ \lbrack D_t t\rbrack^n = 1, \tag{6.33}\\ \lbrack D_{2t} t\rbrack^n = 1, \tag{6.34}\\ \lbrack D_{t}D_t t\rbrack^n = 0 \tag{6.35} \tp \end{align}

The next formulas concern the action of difference operators on a $$t^2$$ term. \begin{align} \lbrack D_t^+ t^2\rbrack^n = (2n+1)\Delta t, \tag{6.36}\\ \lbrack D_t^- t^2\rbrack^n = (2n-1)\Delta t, \tag{6.37}\\ \lbrack D_t t^2\rbrack^n = 2n\Delta t, \tag{6.38}\\ \lbrack D_{2t} t^2\rbrack^n = 2n\Delta t, \tag{6.39}\\ \lbrack D_{t}D_t t^2\rbrack^n = 2, \tag{6.40} \end{align}

Finally, we present formulas for a $$t^3$$ term: \begin{align} \lbrack D_t^+ t^3\rbrack^n &= 3(n\Delta t)^2 + 3n\Delta t^2 + \Delta t^2, \tag{6.41}\\ \lbrack D_t^- t^3\rbrack^n &= 3(n\Delta t)^2 - 3n\Delta t^2 + \Delta t^2, \tag{6.42}\\ \lbrack D_t t^3\rbrack^n &= 3(n\Delta t)^2 + \frac{1}{4}\Delta t^2, \tag{6.43}\\ \lbrack D_{2t} t^3\rbrack^n &= 3(n\Delta t)^2 + \Delta t^2, \tag{6.44}\\ \lbrack D_{t}D_t t^3\rbrack^n &= 6n\Delta t, \tag{6.45} \end{align}

## Software

Application of finite difference operators to polynomials and exponential functions, resulting in the formulas above, can easily be computed by some sympy code (from the file lib.py):

from sympy import *
t, dt, n, w = symbols('t dt n w', real=True)

# Finite difference operators

def D_t_forward(u):
return (u(t + dt) - u(t))/dt

def D_t_backward(u):
return (u(t) - u(t-dt))/dt

def D_t_centered(u):
return (u(t + dt/2) - u(t-dt/2))/dt

def D_2t_centered(u):
return (u(t + dt) - u(t-dt))/(2*dt)

def D_t_D_t(u):
return (u(t + dt) - 2*u(t) + u(t-dt))/(dt**2)

op_list = [D_t_forward, D_t_backward,
D_t_centered, D_2t_centered, D_t_D_t]

def ft1(t):
return t

def ft2(t):
return t**2

def ft3(t):
return t**3

def f_expiwt(t):
return exp(I*w*t)

def f_expwt(t):
return exp(w*t)

func_list = [ft1, ft2, ft3, f_expiwt, f_expwt]


To see the results, one can now make a simple loop over the different type of functions and the various operators associated with them:

for func in func_list:
for op in op_list:
f = func
e = op(f)
e = simplify(expand(e))
print e
if func in [f_expiwt, f_expwt]:
e = e/f(t)
e = e.subs(t, n*dt)
print expand(e)
print factor(simplify(expand(e)))