$$
\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}}
$$
$$
\begin{align}
u'(t_n) &\approx
\lbrack D_tu\rbrack^n = \frac{u^{n+\half} - u^{n-\half}}{\Delta t} \tag{601}\\
u'(t_n) &\approx
\lbrack D_{2t}u\rbrack^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t}
\tag{602}\\
u'(t_n) &=
\lbrack D_t^-u\rbrack^n = \frac{u^{n} - u^{n-1}}{\Delta t}
\tag{603}\\
u'(t_n) &\approx
\lbrack D_t^+u\rbrack^n = \frac{u^{n+1} - u^{n}}{\Delta t}
\tag{604}\\
u'(t_{n+\theta}) &=
\lbrack \bar D_tu\rbrack^{n+\theta} = \frac{u^{n+1} - u^{n}}{\Delta t}
\tag{605}\\
u'(t_n) &\approx
\lbrack D_t^{2-}u\rbrack^n = \frac{3u^{n} - 4u^{n-1} + u^{n-2}}{2\Delta t}
\tag{606}\\
u''(t_n) &\approx
\lbrack D_tD_t u\rbrack^n = \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2}
\tag{607}\\
u(t_{n+\half}) &\approx \lbrack \overline{u}^{t}\rbrack^{n+\half}
= \half(u^{n+1} + u^n)
\tag{608}\\
u(t_{n+\half})^2 &\approx \lbrack \overline{u^2}^{t,g}\rbrack^{n+\half}
= u^{n+1}u^n
\tag{609}\\
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{610}\\
u(t_{n+\theta}) &\approx \lbrack \overline{u}^{t,\theta}\rbrack^{n+\theta}
= \theta u^{n+1} + (1-\theta)u^n ,
\tag{611}\\
&\qquad t_{n+\theta}=\theta t_{n+1} + (1-\theta)t_{n-1}
\tag{612}
\end{align}
$$
Some may wonder why \( \theta \) is absent on the right-hand side
of (605). The fraction is an approximation to the
derivative at the point \( t_{n+\theta}=\theta t_{n+1} + (1-theta) t_{n} \).
$$
\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{613}\\
\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{614}\\
\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{615}\\
\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{616}\\
\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{617}\\
\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{618}\\
\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{619}
\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{620}
\end{align}
$$
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{621}\\
[D_t^+ u]^n &= u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
\tag{622}\\
[D_t^- u]^n &= u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
\tag{623}\\
[D_t u]^n &= u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
\tag{624}\\
[D_{2t} u]^n &= u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}
\tag{625}
\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{626}\\
[D_t^+ u]^n &= u^n\frac{1}{\Delta t}(\exp{(i\omega\Delta t)} - 1),
\tag{627}\\
[D_t^- u]^n &= u^n\frac{1}{\Delta t}(1 - \exp{(-i\omega\Delta t)}),
\tag{628}\\
[D_t u]^n &= u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)},
\tag{629}\\
[D_{2t} u]^n &= u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}
\tag{630}
\tp
\end{align}
$$
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{631}\\
\lbrack D_t^- t\rbrack^n = 1,
\tag{632}\\
\lbrack D_t t\rbrack^n = 1,
\tag{633}\\
\lbrack D_{2t} t\rbrack^n = 1,
\tag{634}\\
\lbrack D_{t}D_t t\rbrack^n = 0
\tag{635}
\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{636}\\
\lbrack D_t^- t^2\rbrack^n = (2n-1)\Delta t,
\tag{637}\\
\lbrack D_t t^2\rbrack^n = 2n\Delta t,
\tag{638}\\
\lbrack D_{2t} t^2\rbrack^n = 2n\Delta t,
\tag{639}\\
\lbrack D_{t}D_t t^2\rbrack^n = 2,
\tag{640}
\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{641}\\
\lbrack D_t^- t^3\rbrack^n &= 3(n\Delta t)^2 - 3n\Delta t^2 + \Delta t^2,
\tag{642}\\
\lbrack D_t t^3\rbrack^n &= 3(n\Delta t)^2 + \frac{1}{4}\Delta t^2,
\tag{643}\\
\lbrack D_{2t} t^3\rbrack^n &= 3(n\Delta t)^2 + \Delta t^2,
\tag{644}\\
\lbrack D_{t}D_t t^3\rbrack^n &= 6n\Delta t,
\tag{645}
\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)))