Can use sympy
and the tools in
truncation_errors.py:
def decay():
u, a = sm.symbols('u a')
diffop = DiffOp(u, independent_variable='t',
num_terms_Taylor_series=3)
D1u = diffop.D(1) # symbol for du/dt
ODE = D1u + a*u # define ODE
# Define schemes
FE = diffop['Dtp'] + a*u
CN = diffop['Dt' ] + a*u
BE = diffop['Dtm'] + a*u
# Residuals (truncation errors)
R = {'FE': FE-ODE, 'BE': BE-ODE, 'CN': CN-ODE}
return R
The returned dictionary becomes
decay: {
'BE': D2u*dt/2 + D3u*dt**2/6,
'FE': -D2u*dt/2 + D3u*dt**2/6,
'CN': D3u*dt**2/24,
}
\( \theta \)-rule: see truncation_errors.py
(long expression, very advantageous
to automate the math!)