Sympy code for FDM expansion of LBM diffusive equation up to Fourth-order
import sympy as sp
# -----------------------------------------------------------------------------
# 1) Symbols (parameters and step sizes)
# -----------------------------------------------------------------------------
dx, dt = sp.symbols('Delta_x Delta_t', positive=True, real=True)
s1, s2, w0, nu = sp.symbols('s_{1} s_{2} w_{0} \\nu', real=True)
alpha1 = sp.simplify(1 - s1/2 - w0*s2/2)
alpha2 = sp.simplify((w0-1)*s2 + 1)
beta1 = sp.simplify(w0*s1*s2/2 - s1*s2/2 - w0*s2/2 + s1/2 +s2 -1)
beta2 = sp.simplify(-w0*s1*s2 + w0*s2 + s1 - 1)
gamma = sp.simplify((s1-1)*(s2-1))
# -----------------------------------------------------------------------------
# 2) Continuous variables and field
# -----------------------------------------------------------------------------
x, t = sp.symbols('x t', real=True)
phi = sp.Function('phi')
# Shorthand for the base (continuous) field value at (x,t)
PHI = phi(x, t)
def taylor_phi(sx=0, st=0, order_x=4, order_t=4):
expr = sp.S(0)
# Keep all mixed terms up to the specified individual orders.
for m in range(order_x + 1):
for n in range(order_t + 1):
# skip the (0,0) term? no, keep it.
term = ( (sx*dx)**m * (st*dt)**n
/ (sp.factorial(m)*sp.factorial(n))
* sp.diff(PHI, x, m, t, n) )
expr += term
return sp.expand(expr)
# -----------------------------------------------------------------------------
# 4) Discrete stencil terms (expanded)
# -----------------------------------------------------------------------------
# phi_x^{t+1}, phi_x^{t}, phi_x^{t-1}, phi_x^{t-2}
phi_t_p1 = taylor_phi(sx=0, st=+1, order_x=0, order_t=4)
display(phi_t_p1)
phi_t_0 = taylor_phi(sx=0, st=0, order_x=0, order_t=0) # exactly phi(x,t)
display(phi_t_0)
phi_t_m1 = taylor_phi(sx=0, st=-1, order_x=0, order_t=4)
display(phi_t_m1)
phi_t_m2 = taylor_phi(sx=0, st=-2, order_x=0, order_t=4)
display(phi_t_m2)
# phi_{x±1}^{t}
phi_xm1_t0 = taylor_phi(sx=-1, st=0, order_x=4, order_t=0)
display(phi_xm1_t0)
phi_xp1_t0 = taylor_phi(sx=+1, st=0, order_x=4, order_t=0)
display(phi_xp1_t0)
# phi_{x±1}^{t-1}
phi_xm1_tm1 = taylor_phi(sx=-1, st=-1, order_x=4, order_t=4)
phi_xp1_tm1 = taylor_phi(sx=+1, st=-1, order_x=4, order_t=4)
\[\displaystyle \frac{\Delta_{t}^{4} \frac{\partial^{4}}{\partial t^{4}} \phi{\left(x,t \right)}}{24} + \frac{\Delta_{t}^{3} \frac{\partial^{3}}{\partial t^{3}} \phi{\left(x,t \right)}}{6} + \frac{\Delta_{t}^{2} \frac{\partial^{2}}{\partial t^{2}} \phi{\left(x,t \right)}}{2} + \Delta_{t} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \phi{\left(x,t \right)}\]
\[\displaystyle \phi{\left(x,t \right)}\]
\[\displaystyle \frac{\Delta_{t}^{4} \frac{\partial^{4}}{\partial t^{4}} \phi{\left(x,t \right)}}{24} - \frac{\Delta_{t}^{3} \frac{\partial^{3}}{\partial t^{3}} \phi{\left(x,t \right)}}{6} + \frac{\Delta_{t}^{2} \frac{\partial^{2}}{\partial t^{2}} \phi{\left(x,t \right)}}{2} - \Delta_{t} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \phi{\left(x,t \right)}\]
\[\displaystyle \frac{2 \Delta_{t}^{4} \frac{\partial^{4}}{\partial t^{4}} \phi{\left(x,t \right)}}{3} - \frac{4 \Delta_{t}^{3} \frac{\partial^{3}}{\partial t^{3}} \phi{\left(x,t \right)}}{3} + 2 \Delta_{t}^{2} \frac{\partial^{2}}{\partial t^{2}} \phi{\left(x,t \right)} - 2 \Delta_{t} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \phi{\left(x,t \right)}\]
\[\displaystyle \frac{\Delta_{x}^{4} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{24} - \frac{\Delta_{x}^{3} \frac{\partial^{3}}{\partial x^{3}} \phi{\left(x,t \right)}}{6} + \frac{\Delta_{x}^{2} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}}{2} - \Delta_{x} \frac{\partial}{\partial x} \phi{\left(x,t \right)} + \phi{\left(x,t \right)}\]
\[\displaystyle \frac{\Delta_{x}^{4} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{24} + \frac{\Delta_{x}^{3} \frac{\partial^{3}}{\partial x^{3}} \phi{\left(x,t \right)}}{6} + \frac{\Delta_{x}^{2} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}}{2} + \Delta_{x} \frac{\partial}{\partial x} \phi{\left(x,t \right)} + \phi{\left(x,t \right)}\]
# -----------------------------------------------------------------------------
# 5) Build the discrete equation (expanded) and move everything to one side
# -----------------------------------------------------------------------------
lhs = phi_t_p1
rhs = (alpha1*phi_xm1_t0 + alpha2*phi_t_0 + alpha1*phi_xp1_t0
+ beta1*phi_xm1_tm1 + beta2*phi_t_m1 + beta1*phi_xp1_tm1
+ gamma * phi_t_m2)
residual = sp.simplify(sp.expand(lhs - rhs)) # = 0 is the modified equation
max_dx=4
max_dt=2
max_total=4
expr = sp.expand(residual)
poly = sp.Poly(expr, dx, dt, domain='EX') # keep symbolic coeffs
kept = sp.S(0)
for (px, pt), coeff in poly.terms():
if (px <= max_dx) and (pt <= max_dt) and ((px + pt) <= max_total):
kept += coeff * dx**px * dt**pt
# display(sp.simplify(residual))
display(sp.simplify(kept.subs(dt**2*dx**2,0)))
\[\displaystyle \frac{\Delta_{t}^{2} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right) \frac{\partial^{2}}{\partial t^{2}} \phi{\left(x,t \right)}}{2} + \frac{\Delta_{t} \Delta_{x}^{2} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right) \frac{\partial^{3}}{\partial x^{2}\partial t} \phi{\left(x,t \right)}}{2} + \Delta_{t} s_{1} s_{2} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \frac{\Delta_{x}^{4} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{24} + \frac{\Delta_{x}^{2} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right) \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}}{2}\]
rel1=sp.simplify(kept.subs(dt**2*dx**2, 0).subs(PHI.diff(t, 2), (1-w0)*((2-s1)/(2*s1))*dx*dx/dt*nu*PHI.diff(x, 4)).subs(PHI.diff(t, x, 2), nu*PHI.diff(x, 4)))
rel1
\[\displaystyle \frac{6 \Delta_{t} \Delta_{x}^{2} \nu s_{1} \left(2 - 3 s_{2}\right) \left(s_{1} - 2\right) \left(w_{0} - 1\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} + 12 \Delta_{t} \Delta_{x}^{2} \nu s_{2} \left(s_{1} - 2\right) \left(w_{0} - 1\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} + s_{1} \left(12 \Delta_{t} \Delta_{x}^{2} \nu \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} + 24 \Delta_{t} s_{1} s_{2} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \Delta_{x}^{4} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} + 12 \Delta_{x}^{2} s_{2} \left(- s_{1} w_{0} + s_{1} + 2 w_{0} - 2\right) \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}\right)}{24 s_{1}}\]
rel2=(sp.collect(sp.simplify(sp.expand(rel1)/(dt*s1*s2)),[PHI.diff(x, 4),PHI.diff(x, 2)])
.subs(sp.expand(dx*dx*(-w0/2+sp.Rational(1,2)+w0/s1-1/s1)/dt),sp.factor(dx*dx*(-w0/2+sp.Rational(1,2)+w0/s1-1/s1)/dt))
.subs(sp.factor(-dx*dx*(-w0/2+sp.Rational(1,2)+w0/s1-1/s1)/dt), nu)
.subs(sp.expand(dx*dx*dx*dx*(-w0/2+sp.Rational(1,2)+w0/s1-1/s1)/(12*dt)), -dx*dx*nu/12)
.subs(sp.expand(dx*dx*nu*(-w0/4 + sp.Rational(1,6) + w0/(2*s2) + 3*w0/(2*s1) - 1/s1 - w0/(s1*s2) - w0/s1**2 + 1/s1**2 ) ), sp.factor(sp.expand(dx*dx*nu*(-w0/4 + sp.Rational(1,6) + w0/(2*s2) + 3*w0/(2*s1) - 1/s1 - w0/(s1*s2) - w0/s1**2 + 1/s1**2 ) ))) )
rel2
\[\displaystyle - \frac{\Delta_{x}^{2} \nu \left(3 s_{1}^{2} s_{2} w_{0} - 2 s_{1}^{2} s_{2} - 6 s_{1}^{2} w_{0} - 18 s_{1} s_{2} w_{0} + 12 s_{1} s_{2} + 12 s_{1} w_{0} + 12 s_{2} w_{0} - 12 s_{2}\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{12 s_{1}^{2} s_{2}} - \nu \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)} + \frac{\partial}{\partial t} \phi{\left(x,t \right)}\]
# from sympy import latex
# print(latex(sp.simplify(kept.subs(dt**2*dx**2,0))))