power series substituting continuity equation

PDE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "from sympy import IndexedBase, symbols, Function, powsimp, Lambda, init_printing\n",
    "init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, j, k, x, r, t = symbols('i j k x r t') #xi, rj, tk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = IndexedBase('b')#rho,q\n",
    "d = IndexedBase('d')#vr,v\n",
    "e = IndexedBase('e')#vx,u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = Function('q')#rho\n",
    "u = Function('u')#vx\n",
    "v = Function('v')#vr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABC8AAAAsBAMAAAC3RI+zAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARHarEIm7zVTvMt2ZImY9RQ3UAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJ8klEQVR4Ad1bb4gkRxV/s9szN7vz55b75Ebl2s0FVA5uOcV8Cjfu6YVDwg6iHJjkHBEShDO36ocVRV0k+aKEWw+jGFT2m/rFW4zBiHc6RPQIiK5IRIkkSwhR0JiJnElQUd971T1V3dOz9ap6uiIWbM/rqvq93/u9quk/U1sAjuXZ9+85Ity7d2+52x3kjAihpBxHmDw4J64QUN+LPlTYMMvKu+HM0iz9FfoKoaQkR5A8FCbHvfJ5gO+4o9wQ3QEsbrtBPHqHUFKOI0wePFJXBHkC4PB+UcMM69o4K/49Q3/FrkIoKccRJg/F2fGpbQ19UI6Yq479vbqHUFKSI0gevJKXB3136/AgXzfr88U/wcuz9jnpL4SSchxh8jCZGZ+abz335JWhD9ABs/jZOz7wL4f+fl1DKCnHESYPftnLo2oPA1yM87UzPl/egcbfZuxzwl0IJSU5guRhIjF+Fb+PAX5S8atk5zJA6yW/+OSoEErKcYTJgzxjB/Y81QP42oE9yjfO3wBYqPx1NYSSchxh8lB+vNjDKwBR1Vf5wxsA6/GMAp7qJoSSchxh8jA1QU4NET4TzuEXutKyPgQ40auUAqd3ACUlOYLkYUZppqtFqz8jZ9Pc0NXi0WmNs6oPoaQkR5A8zCqfFwDeMStf0/w0hjA/gKpLCCXlOBpB8iDKc2d/SrfmXtLw1P1/GGb7PJ091WcPajNjWUk6j93xiwwCnDnASgITStxJrIhyHPk82DUZWbPGRn1FnajjGToUltvS2gdeSC31GQ2y5/psPp1MuootO0nzfBbqzuGhxJ1EgMhlS4Awk5XLgz1xGixiEnVin5/Rnk2rtQe3m+emfSg7huOmaxBtjE8yhjuJOweEIHEPyx1hZs5Fk4hJ1Iki6A7NOLR9CaAe69OM9ZfM2fikOQJ4L58143ElGx4kVo7PZym8lFhJ8kLAishHZUdMCDFc2BNnoEWxiTpRBHhlKCzX8e1uu7AF4FpxfRsvF2e5KZ9PDxIrh5ERFU4VJHkhdukTqXEXYriwazLSYGUix6JO1PEpOhQUenf8WEE9VnV2i+vpnfbQErXl8+lOYucwMqLCqYIkL8QeVj41dsSEEMOFXZNG25nQsaVT7eZ3rsTQOX/sTXAce6/d/OtbYg7nyNuPrJBR+/nLjwPcxXXGYe1xiL4Mc0O8mqzcdkTdNpQv7PSGixdiqGNbMjE8SaLjMfxYwMEZ4XgJcbZHSnRUYFECTz67tuQiRBoWyU+LjEMLSXH8KRwdRkvUiATUG59sbMPy3vwI7sN09rvfa+FtAEv/3q2PsEHPC/BHNo3DcBnmRtDYB3g3/Dn+AbcoX2TiYwl02Q9/0TxJ6rVteETAwRnheNuIOAGkREdlU1L7BpwSkGgh0rA4J+og5NBCDCxIR4fREiaRgG+3Vms70TdhfgMew4fMXv1GbY+iWtx6BG7i8Oh5Ab7Ipj509q9DYwCtLYD3wQd7b+QW9sUWPpZAc5dMzqcnye/au/CAgIMyouJ9GoO9B0iJjsqiBM7uw60CEi1EGhbJT4qQQwtJgfQpHR2eGBImkYDe+hJAewQLffgSxkCjzSWK0pVv/g38Cv6v53+MEsFlWIhhAcE9Ggou7IsteixZHEF9c/PTX93cXPUjgR4yvM7GcW5z84ebm59Q8RLiMpASHZVFCfwK4CEbSUaIKCwwkjUScWSEGHDZ6CRoEZNMwCUcw4UBPLNE6cThj+mIhe8gZNDMAZwY2VK7AVd2OJ8AuIKiCvmiwpMKJwYW/qLxrcWD5GgP1OQ7mIO/Kireo73Ff/DEMBAWJV+BzktOQkAYFslPipDDEJIi8VOaOEaLmEQC6Kq/HsOLfAHGGbCThNTeSIyjeB2B1ycn44/6LrwF33C3sIKGQhXyRaU2wm9sbZdMNTE8Se7ERx0BB2dExXsnzG0rJTqqg5XgUld3JCDRQkAYFslXRcphCEmh+CkdHULLmEQCHiLqLbiXH9mavRMw18ManKd9/gA4B18oePhsrMI93R49fN6PD6HHVF/yRaW+0dwxHz7xWu1F8jZAP3YOzqeK98N4L3yOHj6NqA5WEr0KjQ0BiZoYrFAYFqeCD1IOQ4gGi0eHJ4ZIjURAh77uh+Lu3/mV9OjWRVxaifA/TOCZHZx/ZPw0igF+hkamtFdrX6/TK2nt1caosw/reIFhX2S0V9fwwDOL8+lLcg5fLgQcnE+KF4OF5WEfX66NqFgJRTVFycfhxb6ARAvBr4ooLIonKUIOQ0iKxE/p6DBaxCQRwBf86OSb8VETbxenV9Z+uQ/wo15y8yDjt3/F4D6FfwDGVkyE/OatdLuITh47jj1a38c7yC72IaNz3z7q2cEzdSvxJZk/j37sHJwRjB7L6ZUHT+6gEiMqvg1SVFOUnF65tCQg4YnBQkAYlrH1VMhhCMHX/XQHr3R0GC1iEgqghOJtFlp7ZFFp9vhDGxENORRsxbyWdMRrd2qlxnu4gvOZNHmSWDk4Iyk9fWolaTDj8IqU4POV8fvwBGRSCHNZw8psPRVxZIQYcFHiIEGLmEiBVQB1wp8xQC/T1KmKSmrUYzp7fnLjql6I6VMPKqmhfg2NelzJB08SK8e7NIWytJI0mHFURUpuRZSVxBTCLDZEd2BuwRVxmEISOI+eKHGQoEVMpMAmgFXyA9pH2cTDWt64nSuemNy42thLukZxzihYdvckceJQUaRKJqJK/oEgq4SeKN1JbIj2trkF15lDwZsjEiRKnFJOP8oI1dgEkMP6567i8QyZVPb5aBhpoqE1TJvUZzRIzrtpfWqMadMGbxIXjoQsVZIGA2OjQMm5fw7xOXuQQMc9U2NSiOopQVBWufhxXAV+1RaPDlHJmSQCVPT4clFcusl1oWAr5tR/D3uh2NX0/7qbTuLMEYbEFtbk1lMbIpMyhrf64zrr6Ix7oiFiEnUyvU6zeSumsWw5rV+peiZJF0RLeToAHEIJbz0dLzgfEExxE8N5mbq4/X+oVm3FNJYtq4hNkYyXdqugwDfgh2kLbsVKlmkLrl5wdlXCcF5LcEUG76+2Yuply0oCYJKbxku7lXBACCVq66lecHZUkuxcve4Ie026n+rRxlW9bFlJEEyiF0Qr4YAQStTWU73g7Kgk2bn6qCPsNen+Cq/RGMuWVUShSNRbWhX+2WcIJcnW00ueIhR8/L8Pnl6CwNKtmHrZsgLahITf0ipwr1wGUbI+5C24vvcCBedl6sryMCPHyVZMY9lyRo5NNwnJeEHUbJuZHUTJesxbcOnnJp+i4LxM7QMPirlAG1eTZcvKiJmEF0QrowAIoaTBW095XdlHiYLzMrUPPCiGt2Imy5aVEav9nmqptFqSipWoradqOdZDiILzMrUHOjAktxWzGvb/G5Lc1lPXbJWEu9Ll+/8Xe7jacol3EwwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$r q{\\left (x,r,t \\right )} \\frac{\\partial}{\\partial x} u{\\left (x,r,t \\right )} + r q{\\left (x,r,t \\right )} \\frac{\\partial}{\\partial r} v{\\left (x,r,t \\right )} + r u{\\left (x,r,t \\right )} \\frac{\\partial}{\\partial x} q{\\left (x,r,t \\right )} + r v{\\left (x,r,t \\right )} \\frac{\\partial}{\\partial r} q{\\left (x,r,t \\right )} + r \\frac{\\partial}{\\partial t} q{\\left (x,r,t \\right )} + q{\\left (x,r,t \\right )} v{\\left (x,r,t \\right )}$$"
      ],
      "text/plain": [
       "             ∂                             ∂                             ∂    \n",
       "r⋅q(x, r, t)⋅──(u(x, r, t)) + r⋅q(x, r, t)⋅──(v(x, r, t)) + r⋅u(x, r, t)⋅──(q(\n",
       "             ∂x                            ∂r                            ∂x   \n",
       "\n",
       "                         ∂                  ∂                                 \n",
       "x, r, t)) + r⋅v(x, r, t)⋅──(q(x, r, t)) + r⋅──(q(x, r, t)) + q(x, r, t)⋅v(x, r\n",
       "                         ∂r                 ∂t                                \n",
       "\n",
       "    \n",
       ", t)\n",
       "    "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pde1 = r*q(x,r,t).diff(t,1)+r*u(x,r,t)*q(x,r,t).diff(x,1)+r*q(x,r,t)*u(x,r,t).diff(x,1)+q(x,r,t)*v(x,r,t)+ r*v(x,r,t)*q(x,r,t).diff(r,1)+r*q(x,r,t)*v(x,r,t).diff(r,1)\n",
    "pde1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "qterm = b[i, j, k]*x**i*r**j*t**k\n",
    "uterm = d[i, j, k]*x**i*r**j*t**k\n",
    "vterm = e[i, j, k]*x**i*r**j*t**k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwIAAAAaBAMAAAD2//bgAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJm7MquJRO/dIs12VGbfGimAAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJlElEQVRoBa1ab4hcVxU/b+ffzszOzjMaFVmdMaFi2Fg3iR8CrumoIGtVOgWp2Cx2XEjBtthVRNoa2PWLKG3jIiSRGujTImo/2LW0YB3aTnaTKmXB8Q+o1bJjwQ8Khl1sE9TG9Zxz/513973XnU3uh3fPOffc3++ce9+7c/ewALE2dzgEqLaM7aNKKO5/N4w9ZIwA5DXWNLp2Muoue8lM8q96FojoYcWquxRsVvkphSAYwOXg8srg4YjgKfbIXpwMkMShYlRYBqiZsY8cUNJdcBzgrLECe+UjrRsnO7wrQTKzXJ63OEx/s1V3KZis8vuWFIJgEDnYvLJoOCLos0vm4mSBJI6Vo/x/9ECb+vuV8gjcEgWbegBAeAknO7wrQWKy3AgtDtFb5ZqEAs0u6R0QDCbRnYLLiLIXZ6eIxq88X7mi5Snq9Q6chSfmTdxoFV7CicTdN4nJsqYmRKLfPbKcmSPFZCIYTKLSN0uWEWUvThZK8lhtEybeR0NyBwDWw2Lzk207pbZZumgVmYk1Di9YZpyK8lzeMSD9SUc+PDTN4KxiOyAZbA4ir0ye9TBY6WuPrMXJBFGDwb4jLYACLzhayj1ofo1G2GADOwfVf35JnEO9Yj0kL2roFBxdUrL3zB9tehan7tn/AVL+aHAsMxpRfu4P4pQ4B53NVBIHuU1687esibOK7YBksIlSXmKWne4L52CktKyNiYuzExA1fy8U8eApDDTaHATRCZKn6MGBPQJQbEL902TQbQ6+/ycjs1O1adWYUO/HVKEEbbh3HvW/GZthLkd42QJ4rWUGiL7yziW8ozWdaYfSk8ZPZeV2AFkkg90BzsvOMrO39bgg7yi20Zy+OG8MolHxovdNRzDeRPkYFLrdB7vPqR2obQC8AjCz2rdu5DUNpSPYDimntJVeCO0kJ9R6KI+EMLrobKCZ8cADkvMPfcgNIn2uB5BG4hydxBwAV53lGBzv/rzbHajfgXXFIHNQrtPYiVluPsBfnIIRQW6Ax2Xa4qSCOAwj4Xd6R2QU+Cp8BYCvQ1Nko1djZBlKnVJrNb8YkYkaegWbLRbxQU6JK432Lxgn2fPqjC5BgYn0iGWeZPyxqZmKmUP0jbekkxg/2asdKLmDk7Ny38AkSAb1sdN8ykvOkphuBygiqEdZi5MGIgGV/HwId7SMOf/sy4cgv0gq7wAdNeUOfOLl34XT+aUe2bGRV225wzI+yClxpdGe+Cny6lQ3ofCqgVCYkF9EwxnGz/U+PoYKN6L/czGdxPjJXu1Accna8oso8g7UyHgGJAPngFZ8qTEvMYtt5uF2gCKCGbSnL04aiAGL9Y+GEwehdPFta73Rra0NfDlolB4/eHQWvnzT4d63t7bgFzDbwhvFxR+9wF7BbRE6UCMneHqNL1DKYp6V91xI/J7V6mD4+Ou6NhVnLu3/10GKohoVLuAbuXb7yXuA6PfelUIS3LMSGULRK45qs/ZKxJmprGgHKgde6xGLYlBTOAfOjfKqztprjkDUp9DJW0+ucUSwClmLkwiSEmzpMnTug2Lu2dwiE9KRq78BEvBcti3ojJ8rL1vVCsF/w0ZkNSPcF+XEKWCs+J71lLzQgRE45THTyWrbW+Hvve8ZLZnk/qjSMR6iVxyN/osLA5UZZ8XfAHnFWNQ0l1vjZ4E8Hy0qfwOdD/aPKss0dumLkwiSEmy5XelPwqfKbTzbsH3mxyF19A1wmzQC9sWwuFmKhEGLeH+v93TpptpCI9VcglP0G+JbcdDswGmAf+CvXpyZptj2ebg3fMBoPglXqoKHYU/LFKmeQldd7VEc6xPhuxS+ysruQIxFEbjc1kN4XUKakhjtQL5/mn4q8fb4wNfxmb44PgiviB+sRt4HATKGjXkCBrjtdu4ifuLjjBGwDyA3JVQrYkKNDtzMOtdfygg2sknHpG8tdrvPPNjttnG8CRDiD1acmabYFsIRK9N2xknO4ljh1NpfcV0iduvjk5gtx3ubaGF8lVWF3cirYyTbu9wmIbgsIbkkdne3+91u9+kgwLs7tcKJCJ/pi+ODUFzbglXFtrEmvpV4WsgPCrFNM4xax1c9oVWbsEDp29YIAUanwOyqtpMVm/4GVkimH6wYszfldXJSzSfhSlW9bYZ1H+M4vdpDcwxfuXksymhyexJKsd8vUxKjbyB+fmUsjg9CcW0LViH/FipRDvOQH5SKiJ+lDcC5tt3SsqIQ6gN4v/Ti2x2+sF7udG/Fpnag1oTH1V9ZMeYF/NBcy4vrkk/C9Z16xzmzJDmCK1zmiuEr9ziLhjC5/Y8+X9FMIYl3IHZ+ZSyOD0JxbQuWkfMdGIsa/Q58Q3AKsbhca1m1Fv4axmJLrYbq83joq9KNqipRzaXRh+O1yLfiBLUDHwN4Ecgpznw3/N7SwW/GNuC8UX0SrlSVOwA9VcxRhRpd7WGO8Y3gKt5iEzKLsSgCl9sVePtAQpqSGO/AKBLalrE4PgjFtS1YRv7JnXOfhYXKoCLeNsuAwkh7wqkL/Zvgp061UmGQa1LpBpuqKlHNJdcbP13cZkUXXp3KmTtvWMITIoQ48w1Bj2C4la7mNioDrWDlJE7ClaqxZdgTqSKVKtToag9zYLXr1aKHr9Aki8Z3uX0ODoOENCUx3oEnWiYe7DMWxwehuLYFy8jPb239GwqXbM1W4LNYOTZwpr1rEzcK1Q4EB78IXLoxVSWquQSzF374y21WnMOrU93a2sIdOOCqxQrtsUsWFRHOrzjVJ1GVqpVLA12kUoUaXe1hjlwfVs4nZiZZNJ/Lbc8LcUhTEuMd0KecmpSxOD4Ix+UHa5Bdytckqb8jsKqEdxNX1Umw8upoqu8MSSngRKVqmlGoUGOYJceQFMrdQdqSGO/AUGAOxMSF053RIg8Fmu5MpRts9MeMqLkkWAP7S/Lh+AnE87MfAm7VVqp0kaoeOWbHkY2XNiogwRC9Kc05zS5A3IoIo0VOAxjSjqWbx0LIL+JT1FySrQq6chlP9uGagascAlepwmIO6lioiTEPB+x5W8gYkef0RqoFkXFZ47UgJzNj6WamT/WXmb6ouSRbNcLsjclQ6VYDl39YV6rIFYs5qGOhJsacDrKDEQsZI9rBROliQWRc1ngtyJIlLvex5Aj6nwjcSLLVjQ8pMVzNm4Q6HrB4t71+TUH6REPiJ8Z1XZATA6lFVH+pRfHBZGvcZwhNwb3kzXiJCjU+s+czpMqQ4BPtBsSP67ogJwYyzlWlcW8s2eo57VxVcANvwqBwIgKf2fMZUmVIGAw5y3NPjEsi/x9YIkkdXwHO0AAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$2 i r^{2 j + 1} t^{2 k} x^{2 i - 1} b_{i,j,k} d_{i,j,k} + 2 j r^{2 j} t^{2 k} x^{2 i} b_{i,j,k} e_{i,j,k} + k r^{j + 1} t^{k - 1} x^{i} b_{i,j,k} + r^{2 j} t^{2 k} x^{2 i} b_{i,j,k} e_{i,j,k}$$"
      ],
      "text/plain": [
       "     2⋅j + 1  2⋅k  2⋅i - 1                              2⋅j  2⋅k  2⋅i         \n",
       "2⋅i⋅r       ⋅t   ⋅x       ⋅b[i, j, k]⋅d[i, j, k] + 2⋅j⋅r   ⋅t   ⋅x   ⋅b[i, j, \n",
       "\n",
       "                   j + 1  k - 1  i               2⋅j  2⋅k  2⋅i                \n",
       "k]⋅e[i, j, k] + k⋅r     ⋅t     ⋅x ⋅b[i, j, k] + r   ⋅t   ⋅x   ⋅b[i, j, k]⋅e[i,\n",
       "\n",
       "      \n",
       " j, k]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdeterm1 = powsimp(pde1.replace(q, Lambda((x,r,t), qterm)).doit())\n",
    "pdeterm2 = powsimp(pdeterm1.replace(u, Lambda((x,r,t), uterm)).doit())\n",
    "pdeterm3 = powsimp(pdeterm2.replace(v, Lambda((x,r,t), vterm)).doit())\n",
    "pdeterm3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$2 i r^{2 j + 1} t^{2 k} x^{2 i - 1}  d_{i,j,k} + 2 j r^{2 j} t^{2 k} x^{2 i} e_{i,j,k} + k r^{j + 1} t^{k - 1} x^{i}  + r^{2 j} t^{2 k} x^{2 i}  e_{i,j,k}$$\n",
    "$$(0,0,0)\\implies e_{0,0,0}=0 $$\n",
    "$$(1,0,0)\\implies 2 rx d_{1,0,0}+x^2 e_{1,0,0}=0 $$\n",
    "$$(0,1,0)\\implies r^2 (2 e_{0,1,0}+ e_{0,1,0})=0 $$\n",
    "$$(0,0,1)\\implies  r +t^2 e_{0,0,1}=0 $$\n",
    "$$(1,1,0)\\implies 2 r^3 x d_{1,1,0}+2r^2 x^2e_{1,1,0}+r^2 x^2 e_{1,1,0}=0 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}