Sympy code for FDM expansion of LBM diffusive equation up to Sixth-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=6, order_t=6):
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=6)
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=6)
display(phi_t_m1)
phi_t_m2 = taylor_phi(sx=0, st=-2, order_x=0, order_t=6)
display(phi_t_m2)
# phi_{x±1}^{t}
phi_xm1_t0 = taylor_phi(sx=-1, st=0, order_x=6, order_t=0)
display(phi_xm1_t0)
phi_xp1_t0 = taylor_phi(sx=+1, st=0, order_x=6, order_t=0)
display(phi_xp1_t0)
# phi_{x±1}^{t-1}
phi_xm1_tm1 = taylor_phi(sx=-1, st=-1, order_x=6, order_t=6)
phi_xp1_tm1 = taylor_phi(sx=+1, st=-1, order_x=6, order_t=6)
\[\displaystyle \frac{\Delta_{t}^{6} \frac{\partial^{6}}{\partial t^{6}} \phi{\left(x,t \right)}}{720} + \frac{\Delta_{t}^{5} \frac{\partial^{5}}{\partial t^{5}} \phi{\left(x,t \right)}}{120} + \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}^{6} \frac{\partial^{6}}{\partial t^{6}} \phi{\left(x,t \right)}}{720} - \frac{\Delta_{t}^{5} \frac{\partial^{5}}{\partial t^{5}} \phi{\left(x,t \right)}}{120} + \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{4 \Delta_{t}^{6} \frac{\partial^{6}}{\partial t^{6}} \phi{\left(x,t \right)}}{45} - \frac{4 \Delta_{t}^{5} \frac{\partial^{5}}{\partial t^{5}} \phi{\left(x,t \right)}}{15} + \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}^{6} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{720} - \frac{\Delta_{x}^{5} \frac{\partial^{5}}{\partial x^{5}} \phi{\left(x,t \right)}}{120} + \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}^{6} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{720} + \frac{\Delta_{x}^{5} \frac{\partial^{5}}{\partial x^{5}} \phi{\left(x,t \right)}}{120} + \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=6
max_dt=3
max_total=5
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(kept.subs(dt**2*dx**2,0)))
\[\displaystyle \frac{\Delta_{t}^{3} \left(7 s_{1} s_{2} - 6 s_{1} - 6 s_{2} + 6\right) \frac{\partial^{3}}{\partial t^{3}} \phi{\left(x,t \right)}}{6} + \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}^{4} \left(s_{1} s_{2} w_{0} - s_{1} s_{2} + s_{1} - s_{2} w_{0} + 2 s_{2} - 2\right) \frac{\partial^{5}}{\partial x^{4}\partial t} \phi{\left(x,t \right)}}{24} + \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=(kept.subs(dt**2*dx**2, 0)
.subs(PHI.diff(t, 3), ((1-w0)*((2-s1)/(2*s1)))**2*dx**4/dt**2*nu*PHI.diff(x, 6))
.subs(PHI.diff(t, 2), (1-w0)*((2-s1)/(2*s1))*dx*dx/dt*nu*PHI.diff(x, 4))
.subs(PHI.diff(x, 4, t), nu*PHI.diff(x, 6))
.subs(PHI.diff(t, x, 2), nu*PHI.diff(x, 4)) )
rel1
\[\displaystyle \Delta_{t}^{3} \left(\frac{7 \Delta_{x}^{4} \nu s_{2} \left(1 - w_{0}\right)^{2} \left(2 - s_{1}\right)^{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24 \Delta_{t}^{2} s_{1}} - \frac{\Delta_{x}^{4} \nu \left(1 - w_{0}\right)^{2} \left(2 - s_{1}\right)^{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{4 \Delta_{t}^{2} s_{1}} - \frac{\Delta_{x}^{4} \nu s_{2} \left(1 - w_{0}\right)^{2} \left(2 - s_{1}\right)^{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{4 \Delta_{t}^{2} s_{1}^{2}} + \frac{\Delta_{x}^{4} \nu \left(1 - w_{0}\right)^{2} \left(2 - s_{1}\right)^{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{4 \Delta_{t}^{2} s_{1}^{2}}\right) + \Delta_{t}^{2} \left(- \frac{3 \Delta_{x}^{2} \nu s_{2} \left(1 - w_{0}\right) \left(2 - s_{1}\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{4 \Delta_{t}} + \frac{\Delta_{x}^{2} \nu \left(1 - w_{0}\right) \left(2 - s_{1}\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2 \Delta_{t}} + \frac{\Delta_{x}^{2} \nu s_{2} \left(1 - w_{0}\right) \left(2 - s_{1}\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2 \Delta_{t} s_{1}}\right) + \Delta_{t} \Delta_{x}^{4} \left(\frac{\nu s_{1} s_{2} w_{0} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24} - \frac{\nu s_{1} s_{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24} + \frac{\nu s_{1} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24} - \frac{\nu s_{2} w_{0} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24} + \frac{\nu s_{2} \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{12} - \frac{\nu \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{12}\right) + \Delta_{t} \Delta_{x}^{2} \left(\frac{\nu s_{1} s_{2} w_{0} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2} - \frac{\nu s_{1} s_{2} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2} + \frac{\nu s_{1} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2} - \frac{\nu s_{2} w_{0} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{2} + \nu s_{2} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} - \nu \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}\right) + \Delta_{t} s_{1} s_{2} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \Delta_{x}^{4} \left(- \frac{s_{1} s_{2} w_{0} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{24} + \frac{s_{1} s_{2} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{24} + \frac{s_{2} w_{0} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{12} - \frac{s_{2} \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)}}{12}\right) + \Delta_{x}^{2} \left(- \frac{s_{1} s_{2} w_{0} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}}{2} + \frac{s_{1} s_{2} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}}{2} + s_{2} w_{0} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)} - s_{2} \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)}\right)\]
rel2=(sp.collect(sp.expand(rel1),[PHI.diff(x, 6),PHI.diff(x, 4),PHI.diff(x, 2)]) )
rel2
\[\displaystyle \Delta_{t} s_{1} s_{2} \frac{\partial}{\partial t} \phi{\left(x,t \right)} + \left(- \frac{\Delta_{x}^{2} s_{1} s_{2} w_{0}}{2} + \frac{\Delta_{x}^{2} s_{1} s_{2}}{2} + \Delta_{x}^{2} s_{2} w_{0} - \Delta_{x}^{2} s_{2}\right) \frac{\partial^{2}}{\partial x^{2}} \phi{\left(x,t \right)} + \left(- \frac{\Delta_{t} \Delta_{x}^{2} \nu s_{1} s_{2} w_{0}}{4} + \frac{\Delta_{t} \Delta_{x}^{2} \nu s_{1} s_{2}}{4} + \frac{\Delta_{t} \Delta_{x}^{2} \nu s_{1} w_{0}}{2} + \frac{3 \Delta_{t} \Delta_{x}^{2} \nu s_{2} w_{0}}{2} - \Delta_{t} \Delta_{x}^{2} \nu s_{2} - \Delta_{t} \Delta_{x}^{2} \nu w_{0} - \frac{\Delta_{t} \Delta_{x}^{2} \nu s_{2} w_{0}}{s_{1}} + \frac{\Delta_{t} \Delta_{x}^{2} \nu s_{2}}{s_{1}} - \frac{\Delta_{x}^{4} s_{1} s_{2} w_{0}}{24} + \frac{\Delta_{x}^{4} s_{1} s_{2}}{24} + \frac{\Delta_{x}^{4} s_{2} w_{0}}{12} - \frac{\Delta_{x}^{4} s_{2}}{12}\right) \frac{\partial^{4}}{\partial x^{4}} \phi{\left(x,t \right)} + \left(\frac{7 \Delta_{t} \Delta_{x}^{4} \nu s_{1} s_{2} w_{0}^{2}}{24} - \frac{13 \Delta_{t} \Delta_{x}^{4} \nu s_{1} s_{2} w_{0}}{24} + \frac{\Delta_{t} \Delta_{x}^{4} \nu s_{1} s_{2}}{4} - \frac{\Delta_{t} \Delta_{x}^{4} \nu s_{1} w_{0}^{2}}{4} + \frac{\Delta_{t} \Delta_{x}^{4} \nu s_{1} w_{0}}{2} - \frac{5 \Delta_{t} \Delta_{x}^{4} \nu s_{1}}{24} - \frac{17 \Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}^{2}}{12} + \frac{67 \Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}}{24} - \frac{4 \Delta_{t} \Delta_{x}^{4} \nu s_{2}}{3} + \frac{5 \Delta_{t} \Delta_{x}^{4} \nu w_{0}^{2}}{4} - \frac{5 \Delta_{t} \Delta_{x}^{4} \nu w_{0}}{2} + \frac{7 \Delta_{t} \Delta_{x}^{4} \nu}{6} + \frac{13 \Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}^{2}}{6 s_{1}} - \frac{13 \Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}}{3 s_{1}} + \frac{13 \Delta_{t} \Delta_{x}^{4} \nu s_{2}}{6 s_{1}} - \frac{2 \Delta_{t} \Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}} + \frac{4 \Delta_{t} \Delta_{x}^{4} \nu w_{0}}{s_{1}} - \frac{2 \Delta_{t} \Delta_{x}^{4} \nu}{s_{1}} - \frac{\Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}^{2}}{s_{1}^{2}} + \frac{2 \Delta_{t} \Delta_{x}^{4} \nu s_{2} w_{0}}{s_{1}^{2}} - \frac{\Delta_{t} \Delta_{x}^{4} \nu s_{2}}{s_{1}^{2}} + \frac{\Delta_{t} \Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{2}} - \frac{2 \Delta_{t} \Delta_{x}^{4} \nu w_{0}}{s_{1}^{2}} + \frac{\Delta_{t} \Delta_{x}^{4} \nu}{s_{1}^{2}}\right) \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}\]
rel3=(sp.collect(sp.expand(rel2/(dt*s1*s2)),[PHI.diff(x, 6),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 ) ))))
rel3
\[\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)} + \left(\frac{7 \Delta_{x}^{4} \nu w_{0}^{2}}{24} - \frac{13 \Delta_{x}^{4} \nu w_{0}}{24} + \frac{\Delta_{x}^{4} \nu}{4} - \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{4 s_{2}} + \frac{\Delta_{x}^{4} \nu w_{0}}{2 s_{2}} - \frac{5 \Delta_{x}^{4} \nu}{24 s_{2}} - \frac{17 \Delta_{x}^{4} \nu w_{0}^{2}}{12 s_{1}} + \frac{67 \Delta_{x}^{4} \nu w_{0}}{24 s_{1}} - \frac{4 \Delta_{x}^{4} \nu}{3 s_{1}} + \frac{5 \Delta_{x}^{4} \nu w_{0}^{2}}{4 s_{1} s_{2}} - \frac{5 \Delta_{x}^{4} \nu w_{0}}{2 s_{1} s_{2}} + \frac{7 \Delta_{x}^{4} \nu}{6 s_{1} s_{2}} + \frac{13 \Delta_{x}^{4} \nu w_{0}^{2}}{6 s_{1}^{2}} - \frac{13 \Delta_{x}^{4} \nu w_{0}}{3 s_{1}^{2}} + \frac{13 \Delta_{x}^{4} \nu}{6 s_{1}^{2}} - \frac{2 \Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{2} s_{2}} + \frac{4 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{2} s_{2}} - \frac{2 \Delta_{x}^{4} \nu}{s_{1}^{2} s_{2}} - \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{3}} + \frac{2 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{3}} - \frac{\Delta_{x}^{4} \nu}{s_{1}^{3}} + \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{3} s_{2}} - \frac{2 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{3} s_{2}} + \frac{\Delta_{x}^{4} \nu}{s_{1}^{3} s_{2}}\right) \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)} + \frac{\partial}{\partial t} \phi{\left(x,t \right)}\]
Sixth-Order Term
sixrelex=sp.expand(dx**4*nu*(7*w0**2/24 - 13*w0/24 + sp.Rational(1,4) - w0**2/(4*s2) + w0/(2*s2) - 5/(24*s2) - 17*w0**2/(12*s1) + 67*w0/(24*s1) - 4/(3*s1) + 5*w0**2/(4*s1*s2)
- 5*w0/(2*s1*s2) + 7/(6*s1*s2) + 13*w0**2/(6*s1**2) - 13*w0/(3*s1**2) + 13/(6*s1**2) - 2*w0**2/(s1**2*s2) + 4*w0/(s1**2*s2)
- 2/(s1**2*s2) - w0**2/(s1**3) + 2*w0/(s1**3) - 1/(s1**3) + w0**2/(s1**3*s2) - 2*w0/(s1**3*s2) + 1/(s1**3*s2) ) )
sixrelex
\[\displaystyle \frac{7 \Delta_{x}^{4} \nu w_{0}^{2}}{24} - \frac{13 \Delta_{x}^{4} \nu w_{0}}{24} + \frac{\Delta_{x}^{4} \nu}{4} - \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{4 s_{2}} + \frac{\Delta_{x}^{4} \nu w_{0}}{2 s_{2}} - \frac{5 \Delta_{x}^{4} \nu}{24 s_{2}} - \frac{17 \Delta_{x}^{4} \nu w_{0}^{2}}{12 s_{1}} + \frac{67 \Delta_{x}^{4} \nu w_{0}}{24 s_{1}} - \frac{4 \Delta_{x}^{4} \nu}{3 s_{1}} + \frac{5 \Delta_{x}^{4} \nu w_{0}^{2}}{4 s_{1} s_{2}} - \frac{5 \Delta_{x}^{4} \nu w_{0}}{2 s_{1} s_{2}} + \frac{7 \Delta_{x}^{4} \nu}{6 s_{1} s_{2}} + \frac{13 \Delta_{x}^{4} \nu w_{0}^{2}}{6 s_{1}^{2}} - \frac{13 \Delta_{x}^{4} \nu w_{0}}{3 s_{1}^{2}} + \frac{13 \Delta_{x}^{4} \nu}{6 s_{1}^{2}} - \frac{2 \Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{2} s_{2}} + \frac{4 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{2} s_{2}} - \frac{2 \Delta_{x}^{4} \nu}{s_{1}^{2} s_{2}} - \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{3}} + \frac{2 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{3}} - \frac{\Delta_{x}^{4} \nu}{s_{1}^{3}} + \frac{\Delta_{x}^{4} \nu w_{0}^{2}}{s_{1}^{3} s_{2}} - \frac{2 \Delta_{x}^{4} \nu w_{0}}{s_{1}^{3} s_{2}} + \frac{\Delta_{x}^{4} \nu}{s_{1}^{3} s_{2}}\]
sixrels1=sp.factor(sp.expand(dx**4*nu*(7*w0**2/24 - 13*w0/24 + sp.Rational(1,4) - w0**2/(4*s2) + w0/(2*s2) - 5/(24*s2) - 17*w0**2/(12*s1) + 67*w0/(24*s1) - 4/(3*s1) + 5*w0**2/(4*s1*s2)
- 5*w0/(2*s1*s2) + 7/(6*s1*s2) + 13*w0**2/(6*s1**2) - 13*w0/(3*s1**2) + 13/(6*s1**2) - 2*w0**2/(s1**2*s2) + 4*w0/(s1**2*s2)
- 2/(s1**2*s2) - w0**2/(s1**3) + 2*w0/(s1**3) - 1/(s1**3) + w0**2/(s1**3*s2) - 2*w0/(s1**3*s2) + 1/(s1**3*s2) ) ))
sixrel=sixrels1*24*s1**3*s2/(dx**4*nu)
sixrel
\[\displaystyle 7 s_{1}^{3} s_{2} w_{0}^{2} - 13 s_{1}^{3} s_{2} w_{0} + 6 s_{1}^{3} s_{2} - 6 s_{1}^{3} w_{0}^{2} + 12 s_{1}^{3} w_{0} - 5 s_{1}^{3} - 34 s_{1}^{2} s_{2} w_{0}^{2} + 67 s_{1}^{2} s_{2} w_{0} - 32 s_{1}^{2} s_{2} + 30 s_{1}^{2} w_{0}^{2} - 60 s_{1}^{2} w_{0} + 28 s_{1}^{2} + 52 s_{1} s_{2} w_{0}^{2} - 104 s_{1} s_{2} w_{0} + 52 s_{1} s_{2} - 48 s_{1} w_{0}^{2} + 96 s_{1} w_{0} - 48 s_{1} - 24 s_{2} w_{0}^{2} + 48 s_{2} w_{0} - 24 s_{2} + 24 w_{0}^{2} - 48 w_{0} + 24\]
rel3s2=rel3.subs(sixrelex,sixrels1)
rel3s2
\[\displaystyle \frac{\Delta_{x}^{4} \nu \left(7 s_{1}^{3} s_{2} w_{0}^{2} - 13 s_{1}^{3} s_{2} w_{0} + 6 s_{1}^{3} s_{2} - 6 s_{1}^{3} w_{0}^{2} + 12 s_{1}^{3} w_{0} - 5 s_{1}^{3} - 34 s_{1}^{2} s_{2} w_{0}^{2} + 67 s_{1}^{2} s_{2} w_{0} - 32 s_{1}^{2} s_{2} + 30 s_{1}^{2} w_{0}^{2} - 60 s_{1}^{2} w_{0} + 28 s_{1}^{2} + 52 s_{1} s_{2} w_{0}^{2} - 104 s_{1} s_{2} w_{0} + 52 s_{1} s_{2} - 48 s_{1} w_{0}^{2} + 96 s_{1} w_{0} - 48 s_{1} - 24 s_{2} w_{0}^{2} + 48 s_{2} w_{0} - 24 s_{2} + 24 w_{0}^{2} - 48 w_{0} + 24\right) \frac{\partial^{6}}{\partial x^{6}} \phi{\left(x,t \right)}}{24 s_{1}^{3} s_{2}} - \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)}\]
sp.collect(sixrel.subs(w0,sp.Rational(2,3)),s2)
\[\displaystyle \frac{s_{1}^{3}}{3} + \frac{4 s_{1}^{2}}{3} - \frac{16 s_{1}}{3} + s_{2} \left(\frac{4 s_{1}^{3}}{9} - \frac{22 s_{1}^{2}}{9} + \frac{52 s_{1}}{9} - \frac{8}{3}\right) + \frac{8}{3}\]
fac1=sp.factor(s1**3/3 + 4*s1**2/3 - 16*s1/3 + sp.Rational(8,3))
display(fac1)
fac2=sp.factor(4*s1**3/9 - 22*s1**2/9 + 52*s1/9 - sp.Rational(8,3))
display(fac2)
\[\displaystyle \frac{\left(s_{1} - 2\right) \left(s_{1}^{2} + 6 s_{1} - 4\right)}{3}\]
\[\displaystyle \frac{2 \left(2 s_{1}^{3} - 11 s_{1}^{2} + 26 s_{1} - 12\right)}{9}\]
Fouth-Order Term
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 ) )
\[\displaystyle - \frac{\Delta_{x}^{2} \nu w_{0}}{4} + \frac{\Delta_{x}^{2} \nu}{6} + \frac{\Delta_{x}^{2} \nu w_{0}}{2 s_{2}} + \frac{3 \Delta_{x}^{2} \nu w_{0}}{2 s_{1}} - \frac{\Delta_{x}^{2} \nu}{s_{1}} - \frac{\Delta_{x}^{2} \nu w_{0}}{s_{1} s_{2}} - \frac{\Delta_{x}^{2} \nu w_{0}}{s_{1}^{2}} + \frac{\Delta_{x}^{2} \nu}{s_{1}^{2}}\]
termx=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 ) ))
termx
\[\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)}{12 s_{1}^{2} s_{2}}\]
termx2=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 ) ))/(dx*dx*nu)*12*s1*s1*s2
termx2.subs(w0,sp.Rational(2,3))
\[\displaystyle 4 s_{1}^{2} - 8 s_{1} + 4 s_{2}\]
# from sympy import latex
# print(latex(sp.simplify(kept.subs(dt**2*dx**2,0))))