formulas.qmd

## Finite difference operator notation {#sec-form-fdop}

$$
u'(t_n) \approx
\lbrack D_tu\rbrack^n = \frac{u^{n+\half} - u^{n-\half}}{\Delta t}
$$
$$
u'(t_n) \approx
\lbrack D_{2t}u\rbrack^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t}
$$
$$
u'(t_n) =
\lbrack D_t^-u\rbrack^n = \frac{u^{n} - u^{n-1}}{\Delta t}
$$
$$
u'(t_n) \approx
\lbrack D_t^+u\rbrack^n = \frac{u^{n+1} - u^{n}}{\Delta t}
$$
$$
u'(t_{n+\theta}) =
\lbrack \bar D_tu\rbrack^{n+\theta} = \frac{u^{n+1} - u^{n}}{\Delta t}
$$
$$
u'(t_n) \approx
\lbrack D_t^{2-}u\rbrack^n = \frac{3u^{n} - 4u^{n-1} + u^{n-2}}{2\Delta t}
$$
$$
u''(t_n) \approx
\lbrack D_tD_t u\rbrack^n = \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2}
$$
$$
u(t_{n+\half}) \approx \lbrack \overline{u}^{t}\rbrack^{n+\half}
= \half(u^{n+1} + u^n)
$$
$$
u(t_{n+\half})^2 \approx \lbrack \overline{u^2}^{t,g}\rbrack^{n+\half}
= u^{n+1}u^n
$$
$$
u(t_{n+\half}) \approx \lbrack \overline{u}^{t,h}\rbrack^{n+\half}
= \frac{2}{\frac{1}{u^{n+1}} + \frac{1}{u^n}}
$$
$$
\begin{split}
u(t_{n+\theta}) &\approx \lbrack \overline{u}^{t,\theta}\rbrack^{n+\theta}
= \theta u^{n+1} + (1-\theta)u^n ,
\\
&\qquad t_{n+\theta}=\theta t_{n+1} + (1-\theta)t_{n-1}
\end{split}
$$ {#eq-Dop-fd1-center}

Some may wonder why $\theta$ is absent on the right-hand side
of (@eq-Dop-fd1-center). 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 {#sec-form-truncerr}

$$
\begin{split}
\uex'(t_n) &=
[D_t\uex]^n + R^n = \frac{\uex^{n+\half} - \uex^{n-\half}}{\Delta t} +R^n,\\
R^n &= -\frac{1}{24}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
\end{split}
$$
$$
\begin{split}
\uex'(t_n) &=
[D_{2t}\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n-1}}{2\Delta t} +
R^n,\\
R^n &= -\frac{1}{6}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
\end{split}
$$
$$
\begin{split}
\uex'(t_n) &=
[D_t^-\uex]^n +R^n = \frac{\uex^{n} - \uex^{n-1}}{\Delta t}
+R^n,\\
R^n &= -\half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2)
\end{split}
$$
$$
\begin{split}
\uex'(t_n) &=
[D_t^+\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n}}{\Delta t}
+R^n,\\
R^n &= \half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2)
\end{split}
$$
$$
\begin{split}
\uex'(t_{n+\theta}) &=
[\bar D_t\uex]^{n+\theta} +R^{n+\theta} = \frac{\uex^{n+1} - \uex^{n}}{\Delta t}
+R^{n+\theta},\\
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 +
\\
&\quad {\cal O}(\Delta t^3)
\end{split}
$$
$$
\begin{split}
\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,\\
R^n &= \frac{1}{3}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^3)
\end{split}
$$
$$
\begin{split}
\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,\\
R^n &= -\frac{1}{12}\uex''''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
\end{split}
$$ {#eq-form-trunc-fd1-center}

$$
\begin{split}
\uex(t_{n+\theta}) &= [\overline{\uex}^{t,\theta}]^{n+\theta} +R^{n+\theta}
= \theta \uex^{n+1} + (1-\theta)\uex^n +R^{n+\theta},\\
R^{n+\theta} &= -\half\uex''(t_{n+\theta})\Delta t^2\theta (1-\theta) +
{\cal O}(\Delta t^3)\tp
\end{split}
$$ {#eq-form-trunc-avg-theta}

### Complex exponentials {#sec-form-fdexp}

Let $u^n = \exp{(i\omega n\Delta t)} = e^{i\omega t_n}$.

$$
[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),
$$
$$
[D_t^+ u]^n = u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
$$
$$
[D_t^- u]^n = u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
$$ {#eq-form-exp-fd1-bw}
$$
[D_t u]^n = u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
$$ {#eq-form-exp-fd1c-center}
$$
[D_{2t} u]^n = u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}\tp
$$ {#eq-form-exp-fd2-center}

### Real exponentials

Let $u^n = \exp{(\omega n\Delta t)} = e^{\omega t_n}$.

$$
[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),
$$
$$
[D_t^+ u]^n = u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
$$
$$
[D_t^- u]^n = u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
$$
$$
[D_t u]^n = u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
$$
$$
[D_{2t} u]^n = u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}\tp
$$ {#eq-form-rexp-fd2-center}

## Finite difference formulas for powers of t {#sec-form-fdtn}

The following results are useful when checking if a polynomial term in a
solution fulfills the discrete equation for the numerical method.

$$
\lbrack D_t^+ t\rbrack^n = 1,
$$
$$
\lbrack D_t^- t\rbrack^n = 1,
$$
$$
\lbrack D_t t\rbrack^n = 1,
$$
$$
\lbrack D_{2t} t\rbrack^n = 1,
$$
$$
\lbrack D_{t}D_t t\rbrack^n = 0\tp
$$ {#eq-form-Dop-tn-fw}

The next formulas concern the action of difference operators on a $t^2$ term.

$$
\lbrack D_t^+ t^2\rbrack^n = (2n+1)\Delta t,
$$
$$
\lbrack D_t^- t^2\rbrack^n = (2n-1)\Delta t,
$$
$$
\lbrack D_t t^2\rbrack^n = 2n\Delta t,
$$
$$
\lbrack D_{2t} t^2\rbrack^n = 2n\Delta t,
$$
$$
\lbrack D_{t}D_t t^2\rbrack^n = 2,
$$ {#eq-form-Dop-tn2-fw}

Finally, we present formulas for a $t^3$ term:

$$
\lbrack D_t^+ t^3\rbrack^n = 3(n\Delta t)^2 + 3n\Delta t^2 + \Delta t^2,
$$
$$
\lbrack D_t^- t^3\rbrack^n = 3(n\Delta t)^2 - 3n\Delta t^2 + \Delta t^2,
$$
$$
\lbrack D_t t^3\rbrack^n = 3(n\Delta t)^2 + \frac{1}{4}\Delta t^2,
$$
$$
\lbrack D_{2t} t^3\rbrack^n = 3(n\Delta t)^2 + \Delta t^2,
$$
$$
\lbrack D_{t}D_t t^3\rbrack^n = 6n\Delta t,
$$ {#eq-form-Dop-tn3-fw}

## 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 `src/formulas/lib.py` in this repository):

```python
from sympy import *

t, dt, n, w = symbols("t dt n w", real=True)


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
types of functions and the various operators associated with them:

```python
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)))
```