9 Formulas
9.1 Finite difference operator notation
\[ 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} \tag{9.1}\]
Some may wonder why \(\theta\) is absent on the right-hand side of (9.1). The fraction is an approximation to the derivative at the point \(t_{n+\theta}=\theta t_{n+1} + (1-\theta) t_{n}\).
9.2 Truncation errors of finite difference approximations
\[ \begin{split} u_{\text{e}}'(t_n) &= [D_tu_{\text{e}}]^n + R^n = \frac{u_{\text{e}}^{n+\half} - u_{\text{e}}^{n-\half}}{\Delta t} +R^n,\\ R^n &= -\frac{1}{24}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^4} \end{split} \] \[ \begin{split} u_{\text{e}}'(t_n) &= [D_{2t}u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n-1}}{2\Delta t} + R^n,\\ R^n &= -\frac{1}{6}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^4} \end{split} \] \[ \begin{split} u_{\text{e}}'(t_n) &= [D_t^-u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n} - u_{\text{e}}^{n-1}}{\Delta t} +R^n,\\ R^n &= -\half u_{\text{e}}''(t_n)\Delta t + \Oof{\Delta t^2} \end{split} \] \[ \begin{split} u_{\text{e}}'(t_n) &= [D_t^+u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n}}{\Delta t} +R^n,\\ R^n &= \half u_{\text{e}}''(t_n)\Delta t + \Oof{\Delta t^2} \end{split} \] \[ \begin{split} u_{\text{e}}'(t_{n+\theta}) &= [\bar D_tu_{\text{e}}]^{n+\theta} +R^{n+\theta} = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n}}{\Delta t} +R^{n+\theta},\\ R^{n+\theta} &= -\half(1-2\theta)u_{\text{e}}''(t_{n+\theta})\Delta t + \frac{1}{6}((1 - \theta)^3 - \theta^3)u_{\text{e}}'''(t_{n+\theta})\Delta t^2 + \\ &\quad \Oof{\Delta t^3} \end{split} \] \[ \begin{split} u_{\text{e}}'(t_n) &= [D_t^{2-}u_{\text{e}}]^n +R^n = \frac{3u_{\text{e}}^{n} - 4u_{\text{e}}^{n-1} + u_{\text{e}}^{n-2}}{2\Delta t} +R^n,\\ R^n &= \frac{1}{3}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^3} \end{split} \] \[ \begin{split} u_{\text{e}}''(t_n) &= [D_tD_t u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - 2u_{\text{e}}^{n} + u_{\text{e}}^{n-1}}{\Delta t^2} +R^n,\\ R^n &= -\frac{1}{12}u_{\text{e}}''''(t_n)\Delta t^2 + \Oof{\Delta t^4} \end{split} \tag{9.2}\]
\[ \begin{split} u_{\text{e}}(t_{n+\theta}) &= [\overline{u_{\text{e}}}^{t,\theta}]^{n+\theta} +R^{n+\theta} = \theta u_{\text{e}}^{n+1} + (1-\theta)u_{\text{e}}^n +R^{n+\theta},\\ R^{n+\theta} &= -\half u_{\text{e}}''(t_{n+\theta})\Delta t^2\theta (1-\theta) + \Oof{\Delta t^3}\tp \end{split} \tag{9.3}\]
9.2.1 Complex exponentials
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)}), \tag{9.4}\] \[ [D_t u]^n = u^n\frac{2}{\Delta t}i\sin{\left(\frac{\omega\Delta t}{2}\right)}, \tag{9.5}\] \[ [D_{2t} u]^n = u^n\frac{1}{\Delta t}i\sin{(\omega\Delta t)}\tp \tag{9.6}\]
9.2.2 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 \tag{9.7}\]
9.3 Finite difference formulas for powers of t
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 \tag{9.8}\]
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, \tag{9.9}\]
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, \tag{9.10}\]
9.4 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):
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:
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)))