{ "cells": [ { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Sympy code for FDM expansion of LBM diffusive equation up to Fourth-order" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "outputs": [], "source": [ "import sympy as sp\n", "# -----------------------------------------------------------------------------\n", "# 1) Symbols (parameters and step sizes)\n", "# -----------------------------------------------------------------------------\n", "dx, dt = sp.symbols('Delta_x Delta_t', positive=True, real=True)\n", "s1, s2, w0, nu = sp.symbols('s_{1} s_{2} w_{0} \\\\nu', real=True)\n", "\n", "def d(name, nx=0, nt=0):\n", " n = nx + nt\n", " if n == 1:\n", " num = rf\"\\partial{{{name}}}\"\n", " else:\n", " num = rf\"\\partial^{{{n}}}{{{name}}}\"\n", " den = \"\"\n", " if nx:\n", " den += rf\"\\partial{{x^{{{nx}}}}}\" if nx > 1 else r\"\\partial{x}\"\n", " if nt:\n", " den += rf\"\\partial{{t^{{{nt}}}}}\" if nt > 1 else r\"\\partial{t}\"\n", " return sp.Symbol(rf\"\\frac{{{num}}}{{{den}}}\")\n", "# -------------------1st-Order-Terms----------------------\n", "dphit = d(r\"\\phi\", nt=1); dphix = d(r\"\\phi\", nx=1)\n", "# -------------------2nd-Order-Terms----------------------\n", "dphix2 = d(r\"\\phi\", nx=2); dphixt = d(r\"\\phi\", nx=1, nt=1); dphit2 = d(r\"\\phi\", nt=2)\n", "# -------------------3rd-Order-Terms----------------------\n", "dphix3 = d(r\"\\phi\", nx=3); dphix2t = d(r\"\\phi\", nx=2, nt=1); dphixt2 = d(r\"\\phi\", nx=1, nt=2); dphit3 = d(r\"\\phi\", nt=3)\n", "# -------------------4th-Order-Terms----------------------\n", "dphix4 = d(r\"\\phi\", nx=4); dphix3t = d(r\"\\phi\", nx=3, nt=1); dphix2t2 = d(r\"\\phi\", nx=2, nt=2); dphixt3 = d(r\"\\phi\", nx=1, nt=3); dphit4 = d(r\"\\phi\", nt=4)\n", "\n", "alpha1 = sp.simplify(1 - s1/2 - w0*s2/2)\n", "alpha2 = sp.simplify((w0-1)*s2 + 1)\n", "beta1 = sp.simplify(w0*s1*s2/2 - s1*s2/2 - w0*s2/2 + s1/2 +s2 -1)\n", "beta2 = sp.simplify(-w0*s1*s2 + w0*s2 + s1 - 1)\n", "gamma = sp.simplify((s1-1)*(s2-1))\n", "# -----------------------------------------------------------------------------\n", "# 2) Continuous variables and field\n", "# -----------------------------------------------------------------------------\n", "x, t = sp.symbols('x t', real=True)\n", "phi = sp.Function('phi')\n", "# Shorthand for the base (continuous) field value at (x,t)\n", "PHI = phi(x, t)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def taylor_phi(sx=0, st=0, order_x=4, order_t=4):\n", " expr = sp.S(0)\n", " # Keep all mixed terms up to the specified individual orders.\n", " for m in range(order_x + 1):\n", " for n in range(order_t + 1):\n", " # skip the (0,0) term? no, keep it.\n", " term = ( (sx*dx)**m * (st*dt)**n\n", " / (sp.factorial(m)*sp.factorial(n))\n", " * sp.diff(PHI, x, m, t, n) )\n", " expr += term\n", " return sp.expand(expr)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\phi_x^{t+1} = \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_x^{t} = \\phi{\\left(x,t \\right)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_x^{t-1} = \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_x^{t-2} = \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_{x-1}^{t} = \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_{x+1}^{t} = \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_{x-1}^{t-1} = \\frac{\\Delta_{t}^{4} \\Delta_{x}^{4} \\frac{\\partial^{8}}{\\partial x^{4}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{576} - \\frac{\\Delta_{t}^{4} \\Delta_{x}^{3} \\frac{\\partial^{7}}{\\partial x^{3}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{144} + \\frac{\\Delta_{t}^{4} \\Delta_{x}^{2} \\frac{\\partial^{6}}{\\partial x^{2}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{48} - \\frac{\\Delta_{t}^{4} \\Delta_{x} \\frac{\\partial^{5}}{\\partial x\\partial t^{4}} \\phi{\\left(x,t \\right)}}{24} + \\frac{\\Delta_{t}^{4} \\frac{\\partial^{4}}{\\partial t^{4}} \\phi{\\left(x,t \\right)}}{24} - \\frac{\\Delta_{t}^{3} \\Delta_{x}^{4} \\frac{\\partial^{7}}{\\partial x^{4}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{144} + \\frac{\\Delta_{t}^{3} \\Delta_{x}^{3} \\frac{\\partial^{6}}{\\partial x^{3}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{36} - \\frac{\\Delta_{t}^{3} \\Delta_{x}^{2} \\frac{\\partial^{5}}{\\partial x^{2}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{12} + \\frac{\\Delta_{t}^{3} \\Delta_{x} \\frac{\\partial^{4}}{\\partial x\\partial t^{3}} \\phi{\\left(x,t \\right)}}{6} - \\frac{\\Delta_{t}^{3} \\frac{\\partial^{3}}{\\partial t^{3}} \\phi{\\left(x,t \\right)}}{6} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{4} \\frac{\\partial^{6}}{\\partial x^{4}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{48} - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{3} \\frac{\\partial^{5}}{\\partial x^{3}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{12} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} - \\frac{\\Delta_{t}^{2} \\Delta_{x} \\frac{\\partial^{3}}{\\partial x\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} + \\frac{\\Delta_{t}^{2} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} - \\frac{\\Delta_{t} \\Delta_{x}^{4} \\frac{\\partial^{5}}{\\partial x^{4}\\partial t} \\phi{\\left(x,t \\right)}}{24} + \\frac{\\Delta_{t} \\Delta_{x}^{3} \\frac{\\partial^{4}}{\\partial x^{3}\\partial t} \\phi{\\left(x,t \\right)}}{6} - \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} + \\Delta_{t} \\Delta_{x} \\frac{\\partial^{2}}{\\partial x\\partial t} \\phi{\\left(x,t \\right)} - \\Delta_{t} \\frac{\\partial}{\\partial t} \\phi{\\left(x,t \\right)} + \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle \\phi_{x+1}^{t-1} = \\frac{\\Delta_{t}^{4} \\Delta_{x}^{4} \\frac{\\partial^{8}}{\\partial x^{4}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{576} + \\frac{\\Delta_{t}^{4} \\Delta_{x}^{3} \\frac{\\partial^{7}}{\\partial x^{3}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{144} + \\frac{\\Delta_{t}^{4} \\Delta_{x}^{2} \\frac{\\partial^{6}}{\\partial x^{2}\\partial t^{4}} \\phi{\\left(x,t \\right)}}{48} + \\frac{\\Delta_{t}^{4} \\Delta_{x} \\frac{\\partial^{5}}{\\partial x\\partial t^{4}} \\phi{\\left(x,t \\right)}}{24} + \\frac{\\Delta_{t}^{4} \\frac{\\partial^{4}}{\\partial t^{4}} \\phi{\\left(x,t \\right)}}{24} - \\frac{\\Delta_{t}^{3} \\Delta_{x}^{4} \\frac{\\partial^{7}}{\\partial x^{4}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{144} - \\frac{\\Delta_{t}^{3} \\Delta_{x}^{3} \\frac{\\partial^{6}}{\\partial x^{3}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{36} - \\frac{\\Delta_{t}^{3} \\Delta_{x}^{2} \\frac{\\partial^{5}}{\\partial x^{2}\\partial t^{3}} \\phi{\\left(x,t \\right)}}{12} - \\frac{\\Delta_{t}^{3} \\Delta_{x} \\frac{\\partial^{4}}{\\partial x\\partial t^{3}} \\phi{\\left(x,t \\right)}}{6} - \\frac{\\Delta_{t}^{3} \\frac{\\partial^{3}}{\\partial t^{3}} \\phi{\\left(x,t \\right)}}{6} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{4} \\frac{\\partial^{6}}{\\partial x^{4}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{48} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{3} \\frac{\\partial^{5}}{\\partial x^{3}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{12} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} + \\frac{\\Delta_{t}^{2} \\Delta_{x} \\frac{\\partial^{3}}{\\partial x\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} + \\frac{\\Delta_{t}^{2} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} - \\frac{\\Delta_{t} \\Delta_{x}^{4} \\frac{\\partial^{5}}{\\partial x^{4}\\partial t} \\phi{\\left(x,t \\right)}}{24} - \\frac{\\Delta_{t} \\Delta_{x}^{3} \\frac{\\partial^{4}}{\\partial x^{3}\\partial t} \\phi{\\left(x,t \\right)}}{6} - \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} - \\Delta_{t} \\Delta_{x} \\frac{\\partial^{2}}{\\partial x\\partial t} \\phi{\\left(x,t \\right)} - \\Delta_{t} \\frac{\\partial}{\\partial t} \\phi{\\left(x,t \\right)} + \\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)}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from IPython.display import display, Math\n", "# phi_x^{t+1}, phi_x^{t}, phi_x^{t-1}, phi_x^{t-2}\n", "phi_t_p1 = taylor_phi(sx=0, st=+1, order_x=0, order_t=4)\n", "display( Math(r\"\\phi_x^{t+1} = \" + sp.latex(phi_t_p1)))\n", "phi_t_0 = taylor_phi(sx=0, st=0, order_x=0, order_t=0) # exactly phi(x,t)\n", "display( Math(r\"\\phi_x^{t} = \" + sp.latex(phi_t_0)))\n", "phi_t_m1 = taylor_phi(sx=0, st=-1, order_x=0, order_t=4)\n", "display( Math(r\"\\phi_x^{t-1} = \" + sp.latex(phi_t_m1)))\n", "phi_t_m2 = taylor_phi(sx=0, st=-2, order_x=0, order_t=4)\n", "display( Math(r\"\\phi_x^{t-2} = \" + sp.latex(phi_t_m2)))\n", "\n", "# phi_{x±1}^{t}\n", "phi_xm1_t0 = taylor_phi(sx=-1, st=0, order_x=4, order_t=0)\n", "display( Math(r\"\\phi_{x-1}^{t} = \" + sp.latex(phi_xm1_t0)))\n", "phi_xp1_t0 = taylor_phi(sx=+1, st=0, order_x=4, order_t=0)\n", "display( Math(r\"\\phi_{x+1}^{t} = \" + sp.latex(phi_xp1_t0)))\n", "\n", "# phi_{x±1}^{t-1}\n", "phi_xm1_tm1 = taylor_phi(sx=-1, st=-1, order_x=4, order_t=4)\n", "display( Math(r\"\\phi_{x-1}^{t-1} = \" + sp.latex(phi_xm1_tm1)))\n", "phi_xp1_tm1 = taylor_phi(sx=+1, st=-1, order_x=4, order_t=4)\n", "display( Math(r\"\\phi_{x+1}^{t-1} = \" + sp.latex(phi_xp1_tm1)))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# -----------------------------------------------------------------------------\n", "# 5) Build the discrete equation (expanded) and move everything to one side\n", "# -----------------------------------------------------------------------------\n", "lhs = phi_t_p1\n", "\n", "rhs = (alpha1*phi_xm1_t0 + alpha2*phi_t_0 + alpha1*phi_xp1_t0\n", " + beta1*phi_xm1_tm1 + beta2*phi_t_m1 + beta1*phi_xp1_tm1\n", " + gamma * phi_t_m2)\n", "\n", "residual = sp.simplify(sp.expand(lhs - rhs)) # = 0 is the modified equation" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\Delta_{t}^{2} \\Delta_{x}^{2} \\left(- \\frac{s_{1} s_{2} w_{0} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} + \\frac{s_{1} s_{2} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} - \\frac{s_{1} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} + \\frac{s_{2} w_{0} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{4} - \\frac{s_{2} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} + \\frac{\\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2}\\right) + \\Delta_{t}^{2} \\left(- \\frac{3 s_{1} s_{2} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} + s_{1} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)} + s_{2} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}\\right) + \\Delta_{t} \\Delta_{x}^{2} \\left(\\frac{s_{1} s_{2} w_{0} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} - \\frac{s_{1} s_{2} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} + \\frac{s_{1} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} - \\frac{s_{2} w_{0} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)}}{2} + s_{2} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)} - \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\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)$" ], "text/plain": [ "Delta_t**2*Delta_x**2*(-s_{1}*s_{2}*w_{0}*Derivative(phi(x, t), (t, 2), (x, 2))/4 + s_{1}*s_{2}*Derivative(phi(x, t), (t, 2), (x, 2))/4 - s_{1}*Derivative(phi(x, t), (t, 2), (x, 2))/4 + s_{2}*w_{0}*Derivative(phi(x, t), (t, 2), (x, 2))/4 - s_{2}*Derivative(phi(x, t), (t, 2), (x, 2))/2 + Derivative(phi(x, t), (t, 2), (x, 2))/2) + Delta_t**2*(-3*s_{1}*s_{2}*Derivative(phi(x, t), (t, 2))/2 + s_{1}*Derivative(phi(x, t), (t, 2)) + s_{2}*Derivative(phi(x, t), (t, 2))) + Delta_t*Delta_x**2*(s_{1}*s_{2}*w_{0}*Derivative(phi(x, t), t, (x, 2))/2 - s_{1}*s_{2}*Derivative(phi(x, t), t, (x, 2))/2 + s_{1}*Derivative(phi(x, t), t, (x, 2))/2 - s_{2}*w_{0}*Derivative(phi(x, t), t, (x, 2))/2 + s_{2}*Derivative(phi(x, t), t, (x, 2)) - Derivative(phi(x, t), t, (x, 2))) + Delta_t*s_{1}*s_{2}*Derivative(phi(x, t), t) + Delta_x**4*(-s_{1}*s_{2}*w_{0}*Derivative(phi(x, t), (x, 4))/24 + s_{1}*s_{2}*Derivative(phi(x, t), (x, 4))/24 + s_{2}*w_{0}*Derivative(phi(x, t), (x, 4))/12 - s_{2}*Derivative(phi(x, t), (x, 4))/12) + Delta_x**2*(-s_{1}*s_{2}*w_{0}*Derivative(phi(x, t), (x, 2))/2 + s_{1}*s_{2}*Derivative(phi(x, t), (x, 2))/2 + s_{2}*w_{0}*Derivative(phi(x, t), (x, 2)) - s_{2}*Derivative(phi(x, t), (x, 2)))" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "max_dx=4\n", "max_dt=2\n", "max_total=4\n", "expr = sp.expand(residual)\n", "poly = sp.Poly(expr, dx, dt, domain='EX') # keep symbolic coeffs\n", "kept = sp.S(0)\n", "\n", "for (px, pt), coeff in poly.terms():\n", " if (px <= max_dx) and (pt <= max_dt) and ((px + pt) <= max_total):\n", " kept += coeff * dx**px * dt**pt\n", "kept" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} s_{1} s_{2} w_{0}}{4} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} s_{1} s_{2}}{4} - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} s_{1}}{4} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} s_{2} w_{0}}{4} - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} s_{2}}{2} + \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}}{2} - \\frac{3 \\Delta_{t}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}} s_{1} s_{2}}{2} + \\Delta_{t}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}} s_{1} + \\Delta_{t}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}} s_{2} + \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} s_{1} s_{2} w_{0}}{2} - \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} s_{1} s_{2}}{2} + \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} s_{1}}{2} - \\frac{\\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} s_{2} w_{0}}{2} + \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} s_{2} - \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} + \\Delta_{t} \\frac{\\partial{\\phi}}{\\partial{t}} s_{1} s_{2} - \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} s_{1} s_{2} w_{0}}{24} + \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} s_{1} s_{2}}{24} + \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} s_{2} w_{0}}{12} - \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} s_{2}}{12} - \\frac{\\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} s_{1} s_{2} w_{0}}{2} + \\frac{\\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} s_{1} s_{2}}{2} + \\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} s_{2} w_{0} - \\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} s_{2}$" ], "text/plain": [ "-Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*s_{1}*s_{2}*w_{0}/4 + Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*s_{1}*s_{2}/4 - Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*s_{1}/4 + Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*s_{2}*w_{0}/4 - Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*s_{2}/2 + Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}/2 - 3*Delta_t**2*\\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}}*s_{1}*s_{2}/2 + Delta_t**2*\\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}}*s_{1} + Delta_t**2*\\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}}*s_{2} + Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*s_{1}*s_{2}*w_{0}/2 - Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*s_{1}*s_{2}/2 + Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*s_{1}/2 - Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*s_{2}*w_{0}/2 + Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*s_{2} - Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} + Delta_t*\\frac{\\partial{\\phi}}{\\partial{t}}*s_{1}*s_{2} - Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*s_{1}*s_{2}*w_{0}/24 + Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*s_{1}*s_{2}/24 + Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*s_{2}*w_{0}/12 - Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*s_{2}/12 - Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*s_{1}*s_{2}*w_{0}/2 + Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*s_{1}*s_{2}/2 + Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*s_{2}*w_{0} - Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*s_{2}" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "term1 = sp.expand(kept).subs({\n", " # ------------------------------------------------4th order--------------------------------------------------------------------------------------------\n", " PHI.diff(x, 4):dphix4,PHI.diff(x, 3, t):dphix3t,PHI.diff(x, 2, t, 2):dphix2t2,PHI.diff(x, t, 3):dphixt3,PHI.diff(t, 4):dphit4,\n", " # ------------------------------------------------3rd order--------------------------------------------------------------------------------------------\n", " PHI.diff(x, 3):dphix3, PHI.diff(x, 2, t):dphix2t, PHI.diff(x, t, 2):dphixt2, PHI.diff(t, 3):dphit3,\n", " # ------------------------------------------------2nd order--------------------------------------------------------------------------------------------\n", " PHI.diff(x, 2):dphix2, PHI.diff(x, t):dphixt, PHI.diff(t, 2):dphit2,\n", " # ------------------------------------------------1st order--------------------------------------------------------------------------------------------\n", " PHI.diff(x):dphix, PHI.diff(t):dphit,\n", "})\n", "display(term1)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{\\partial\\phi}{\\partial t}= \\Delta_{t} \\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}} \\left(\\frac{3}{2} - \\frac{1}{s_{2}} - \\frac{1}{s_{1}}\\right) + \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} \\left(- \\frac{w_{0}}{2} + \\frac{1}{2} - \\frac{1}{2 s_{2}} + \\frac{w_{0}}{2 s_{1}} - \\frac{1}{s_{1}} + \\frac{1}{s_{1} s_{2}}\\right) + \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} \\left(\\frac{w_{0}}{24} - \\frac{1}{24} - \\frac{w_{0}}{12 s_{1}} + \\frac{1}{12 s_{1}}\\right)}{\\Delta_{t}} + \\frac{\\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} \\left(\\frac{w_{0}}{2} - \\frac{1}{2} - \\frac{w_{0}}{s_{1}} + \\frac{1}{s_{1}}\\right)}{\\Delta_{t}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "term2 = sp.collect(sp.expand(-term1/(dt*s1*s2))+ dphit,[dphit,dphix,\n", " dphix2*dx*dx/dt,dphixt,dphit2*dt,\n", " dphix2t*dx*dx,dphixt2,dphix3,dphit3*dt**2,\n", " dphix2t2*dt*dx**2,dphix4*dx**4/dt,\n", " ]).subs(dphix2t2,0)\n", "display( Math(r\"\\frac{\\partial\\phi}{\\partial t}= \" + sp.latex(term2)))" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{\\partial\\phi}{\\partial t}= \\frac{\\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} \\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)}{24 s_{1}^{2} s_{2}} + \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} \\nu$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nut=dx*dx*(w0/2-sp.Rational(1,2)-w0/s1+1/s1)/dt\n", "term3 = sp.collect(sp.expand(term2.subs(nut,nu)\n", " .subs(dphit2,nut*nut*dphix4)\n", " .subs(dphix2t,nut*dphix4)), [dphix2*dx**2/dt,dphix4*dx**4/dt])\n", "dphix4_term3=sp.simplify(sp.factor(term3.coeff(dphix4)).subs(dt*s1,(w0-1)*(2-s1)*dx**2/nu))\n", "# display(dphix4_term3)\n", "display( Math(r\"\\frac{\\partial\\phi}{\\partial t}= \" + sp.latex(term3.subs(term3.coeff(dphix4),dphix4_term3))))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# from sympy import latex\n", "# print(latex(sp.simplify(kept.subs(dt**2*dx**2,0))))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }