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} \).
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} $$
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} $$
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} $$
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)))