{ "cells": [ { "cell_type": "markdown", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# Sympy code for FDM expansion of LBM diffusive equation" ] }, { "cell_type": "code", "execution_count": 3, "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", "omega, a = sp.symbols('omega a', 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", "Omega = 1 - omega\n", "alpha1 = sp.simplify(Omega + a*omega)\n", "alpha2 = sp.simplify(Omega + (1 - 2*a)*omega)\n", "beta1 = sp.simplify(-Omega*(Omega + (1 - a)*omega))\n", "beta2 = sp.simplify(-Omega*(Omega + 2*a*omega))\n", "\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": 4, "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": 5, "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": 6, "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", " + Omega**2 * phi_t_m2)\n", "\n", "residual = sp.simplify(sp.expand(lhs - rhs)) # = 0 is the modified equation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\Delta_{t}^{2} \\Delta_{x}^{2} \\left(\\frac{a \\omega^{2} \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} - \\frac{a \\omega \\frac{\\partial^{4}}{\\partial x^{2}\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} - \\frac{\\omega \\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 \\omega^{2} \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}}{2} + 2 \\omega \\frac{\\partial^{2}}{\\partial t^{2}} \\phi{\\left(x,t \\right)}\\right) + \\Delta_{t} \\Delta_{x}^{2} \\left(- a \\omega^{2} \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)} + a \\omega \\frac{\\partial^{3}}{\\partial x^{2}\\partial t} \\phi{\\left(x,t \\right)} + \\omega \\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} \\omega^{2} \\frac{\\partial}{\\partial t} \\phi{\\left(x,t \\right)} + \\Delta_{x}^{4} \\left(\\frac{a \\omega^{2} \\frac{\\partial^{4}}{\\partial x^{4}} \\phi{\\left(x,t \\right)}}{12} - \\frac{a \\omega \\frac{\\partial^{4}}{\\partial x^{4}} \\phi{\\left(x,t \\right)}}{6}\\right) + \\Delta_{x}^{2} \\left(a \\omega^{2} \\frac{\\partial^{2}}{\\partial x^{2}} \\phi{\\left(x,t \\right)} - 2 a \\omega \\frac{\\partial^{2}}{\\partial x^{2}} \\phi{\\left(x,t \\right)}\\right)$" ], "text/plain": [ "Delta_t**2*Delta_x**2*(a*omega**2*Derivative(phi(x, t), (t, 2), (x, 2))/2 - a*omega*Derivative(phi(x, t), (t, 2), (x, 2))/2 - omega*Derivative(phi(x, t), (t, 2), (x, 2))/2 + Derivative(phi(x, t), (t, 2), (x, 2))/2) + Delta_t**2*(-3*omega**2*Derivative(phi(x, t), (t, 2))/2 + 2*omega*Derivative(phi(x, t), (t, 2))) + Delta_t*Delta_x**2*(-a*omega**2*Derivative(phi(x, t), t, (x, 2)) + a*omega*Derivative(phi(x, t), t, (x, 2)) + omega*Derivative(phi(x, t), t, (x, 2)) - Derivative(phi(x, t), t, (x, 2))) + Delta_t*omega**2*Derivative(phi(x, t), t) + Delta_x**4*(a*omega**2*Derivative(phi(x, t), (x, 4))/12 - a*omega*Derivative(phi(x, t), (x, 4))/6) + Delta_x**2*(a*omega**2*Derivative(phi(x, t), (x, 2)) - 2*a*omega*Derivative(phi(x, t), (x, 2)))" ] }, "execution_count": 7, "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": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} a \\omega^{2}}{2} - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} a \\omega}{2} - \\frac{\\Delta_{t}^{2} \\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}} \\omega}{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}}} \\omega^{2}}{2} + 2 \\Delta_{t}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}} \\omega - \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} a \\omega^{2} + \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} a \\omega + \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} \\omega - \\Delta_{t} \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} + \\Delta_{t} \\frac{\\partial{\\phi}}{\\partial{t}} \\omega^{2} + \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} a \\omega^{2}}{12} - \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} a \\omega}{6} + \\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} a \\omega^{2} - 2 \\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} a \\omega$" ], "text/plain": [ "Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*a*omega**2/2 - Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*a*omega/2 - Delta_t**2*Delta_x**2*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{2}}\\partial{t^{2}}}*omega/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}}}*omega**2/2 + 2*Delta_t**2*\\frac{\\partial^{2}{\\phi}}{\\partial{t^{2}}}*omega - Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*a*omega**2 + Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*a*omega + Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}}*omega - Delta_t*Delta_x**2*\\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} + Delta_t*\\frac{\\partial{\\phi}}{\\partial{t}}*omega**2 + Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*a*omega**2/12 - Delta_x**4*\\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}}*a*omega/6 + Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*a*omega**2 - 2*Delta_x**2*\\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}}*a*omega" ] }, "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": 7, "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{2}{\\omega}\\right) + \\Delta_{x}^{2} \\frac{\\partial^{3}{\\phi}}{\\partial{x^{2}}\\partial{t}} \\left(a - \\frac{a}{\\omega} - \\frac{1}{\\omega} + \\frac{1}{\\omega^{2}}\\right) + \\frac{\\Delta_{x}^{4} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} \\left(- \\frac{a}{12} + \\frac{a}{6 \\omega}\\right)}{\\Delta_{t}} + \\frac{\\Delta_{x}^{2} \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} \\left(- a + \\frac{2 a}{\\omega}\\right)}{\\Delta_{t}}$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "term2 = sp.collect(sp.expand(-term1/(dt*omega**2))+ 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": 8, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{\\partial\\phi}{\\partial t}= \\frac{\\Delta_{x}^{2} \\frac{\\partial^{4}{\\phi}}{\\partial{x^{4}}} \\nu \\left(- 6 a \\omega^{2} + 48 a \\omega - 48 a + \\omega^{2} - 12 \\omega + 12\\right)}{12 \\omega^{2}} + \\frac{\\partial^{2}{\\phi}}{\\partial{x^{2}}} \\nu$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nus = sp.symbols('\\\\nu')\n", "nut=dx**2/dt*(-a+(2*a/omega))\n", "term3 = sp.collect(sp.expand(term2.subs(nut,nus)\n", " .subs(dphit2,nut*nut*dphix4)\n", " .subs(dphix2t,nut*dphix4)), [dphix4*dx**4/dt])\n", "dphix4_term3=sp.simplify(sp.factor(term3.coeff(dphix4)).subs(dt*omega,a*(2-omega)*dx**2/nus))\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": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }