Let $u^n = \exp{(\omega n\Delta t)} = e^{\omega t_n}$.
The following results are useful when checking if a polynomial term in a solution fulfills the discrete equation for the numerical method.
The next formulas concern the action of difference operators on a $t^2$ term.
Finally, we present formulas for a $t^3$ term:
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 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))))
1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 dt + 2*t 2*dt*n + dt dt*(2*n + 1) -dt + 2*t 2*dt*n - dt dt*(2*n - 1) 2*t 2*dt*n 2*dt*n 2*t 2*dt*n 2*dt*n 2 2 2 dt**2 + 3*dt*t + 3*t**2 3*dt**2*n**2 + 3*dt**2*n + dt**2 dt**2*(3*n**2 + 3*n + 1) dt**2 - 3*dt*t + 3*t**2 3*dt**2*n**2 - 3*dt**2*n + dt**2 dt**2*(3*n**2 - 3*n + 1) dt**2/4 + 3*t**2 3*dt**2*n**2 + dt**2/4 dt**2*(12*n**2 + 1)/4 dt**2 + 3*t**2 3*dt**2*n**2 + dt**2 dt**2*(3*n**2 + 1) 6*t 6*dt*n 6*dt*n (-exp(I*t*w) + exp(I*w*(dt + t)))/dt exp(I*dt*w)/dt - 1/dt (exp(I*dt*w) - 1)/dt (-exp(I*t*w) + exp(I*w*(dt + t)))*exp(-I*dt*w)/dt 1/dt - exp(-I*dt*w)/dt (exp(I*dt*w) - 1)*exp(-I*dt*w)/dt (-exp(I*t*w) + exp(I*w*(dt + t)))*exp(-I*dt*w/2)/dt exp(I*dt*w/2)/dt - exp(-I*dt*w/2)/dt 2*I*sin(dt*w/2)/dt (-exp(I*t*w) + exp(I*w*(2*dt + t)))*exp(-I*dt*w)/(2*dt) exp(I*dt*w)/(2*dt) - exp(-I*dt*w)/(2*dt) I*sin(dt*w)/dt (exp(2*I*dt*w) - 2*exp(I*dt*w) + 1)*exp(I*w*(-dt + t))/dt**2 exp(I*dt*w)/dt**2 - 2/dt**2 + exp(-I*dt*w)/dt**2 (exp(I*dt*w) - 1)**2*exp(-I*dt*w)/dt**2 (-exp(t*w) + exp(w*(dt + t)))/dt exp(dt*w)/dt - 1/dt (exp(dt*w) - 1)/dt (-exp(t*w) + exp(w*(dt + t)))*exp(-dt*w)/dt 1/dt - exp(-dt*w)/dt (exp(dt*w) - 1)*exp(-dt*w)/dt (-exp(t*w) + exp(w*(dt + t)))*exp(-dt*w/2)/dt exp(dt*w/2)/dt - exp(-dt*w/2)/dt 2*sinh(dt*w/2)/dt (-exp(t*w) + exp(w*(2*dt + t)))*exp(-dt*w)/(2*dt) exp(dt*w)/(2*dt) - exp(-dt*w)/(2*dt) sinh(dt*w)/dt (exp(2*dt*w) - 2*exp(dt*w) + 1)*exp(w*(-dt + t))/dt**2 exp(dt*w)/dt**2 - 2/dt**2 + exp(-dt*w)/dt**2 (exp(dt*w) - 1)**2*exp(-dt*w)/dt**2