Poly in s(n) and x_0

Author: hirish99

sumpy/recurrence.py

@@ -401,13 +401,19 @@ def get_shifted_recurrence_exp_from_pde(pde: LinearPDESystemOperator) -> sp.Expr
 
     idx_l, terms = _extract_idx_terms_from_recurrence(r)
 
-    r_ret = r
+    # How much do we need to shift the recurrence relation
+    shift_idx = max(idx_l)
 
     n = sp.symbols("n")
+    r = r.subs(n, n-shift_idx)
+
+    idx_l, terms = _extract_idx_terms_from_recurrence(r)
+
+    r_ret = r
     for i in range(len(idx_l)):
         r_ret = r_ret.subs(terms[i], (-1)**(n+idx_l[i])*terms[i])
 
-    return r_ret
+    return r_ret, (max(idx_l)+1-min(idx_l))
 
 
 def shift_recurrence(r: sp.Expr) -> sp.Expr:

test/taylor_recurrence.ipynb

@@ -4,12 +4,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Part 1: Generalizing the Recurrence:"
+    "# Generalizing a Taylor Recurrence"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -35,7 +35,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 29,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -46,7 +46,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -59,7 +59,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -75,7 +75,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 32,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -93,11 +93,121 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 33,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "4"
+      ]
+     },
+     "execution_count": 33,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "recur, order = get_shifted_recurrence_exp_from_pde(laplace2d)\n",
+    "order"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle 1.55431223447522 \\cdot 10^{-15}$"
+      ],
+      "text/plain": [
+       "1.55431223447522e-15"
+      ]
+     },
+     "execution_count": 34,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "#Sanity check that recurrence is correct\n",
+    "derivs_lap = compute_derivatives(5)\n",
+    "exp = recur.subs(n, 4)\n",
+    "exp.subs(s(4), derivs_lap[4]).subs(s(3), derivs_lap[3]).subs(s(2), derivs_lap[2]).subs(s(1), derivs_lap[1]).subs(var[0],np.random.rand()).subs(var[1],np.random.rand())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 2: We need to arrange the terms in the recurrence as a polynomial in $x_0$, $s(n)$\n",
+    "$$\n",
+    "table[i, j]\n",
+    "$$\n",
+    "Where $i = 0$ represents the coefficient attached to $s(n)$ and $i = 1$ represents $s(n-1)$, etc. and the $j$ is for the polynomial in $x_0$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\operatorname{Poly}{\\left( \\left(\\left(-1\\right)^{n} x_{0}^{3} + \\left(-1\\right)^{n} x_{0} x_{1}^{2}\\right) s{\\left(n \\right)} + \\left(- 3 \\left(-1\\right)^{n} n x_{0}^{2} - \\left(-1\\right)^{n} n x_{1}^{2} + 5 \\left(-1\\right)^{n} x_{0}^{2} + 3 \\left(-1\\right)^{n} x_{1}^{2}\\right) s{\\left(n - 1 \\right)} + \\left(3 \\left(-1\\right)^{n} n^{2} x_{0} - 13 \\left(-1\\right)^{n} n x_{0} + 14 \\left(-1\\right)^{n} x_{0}\\right) s{\\left(n - 2 \\right)} + \\left(- \\left(-1\\right)^{n} n^{3} + 8 \\left(-1\\right)^{n} n^{2} - 21 \\left(-1\\right)^{n} n + 18 \\left(-1\\right)^{n}\\right) s{\\left(n - 3 \\right)}, s{\\left(n \\right)}, s{\\left(n - 1 \\right)}, s{\\left(n - 2 \\right)}, s{\\left(n - 3 \\right)}, domain=\\mathbb{Z}\\left[n, x_{0}, \\left(-1\\right)^{n}, x_{1}\\right] \\right)}$"
+      ],
+      "text/plain": [
+       "Poly(((-1)**n*x0**3 + (-1)**n*x0*x1**2)*(s(n)) + (-3*(-1)**n*n*x0**2 - (-1)**n*n*x1**2 + 5*(-1)**n*x0**2 + 3*(-1)**n*x1**2)*(s(n - 1)) + (3*(-1)**n*n**2*x0 - 13*(-1)**n*n*x0 + 14*(-1)**n*x0)*(s(n - 2)) + (-(-1)**n*n**3 + 8*(-1)**n*n**2 - 21*(-1)**n*n + 18*(-1)**n)*(s(n - 3)), s(n), s(n - 1), s(n - 2), s(n - 3), domain='ZZ[n,x0,(-1)**n,x1]')"
+      ]
+     },
+     "execution_count": 59,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "recur\n",
+    "poly_in_s_n = sp.poly(recur, [s(n-i) for i in range(order)])\n",
+    "poly_in_s_n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "recur = get_shifted_recurrence_exp_from_pde(laplace2d)"
+    "coeff_s_n = [poly_in_s_n.coeff_monomial(poly_in_s_n.gens[i]) for i in range(order)]\n",
+    "\n",
+    "table = []\n",
+    "for i in range(len(coeff_s_n)):\n",
+    "    table.append(sp.poly(coeff_s_n[i], var[0]).all_coeffs()[::-1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 58,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[0, (-1)**n*x1**2, 0, (-1)**n],\n",
+       " [-(-1)**n*n*x1**2 + 3*(-1)**n*x1**2, 0, -3*(-1)**n*n + 5*(-1)**n],\n",
+       " [0, 3*(-1)**n*n**2 - 13*(-1)**n*n + 14*(-1)**n],\n",
+       " [-(-1)**n*n**3 + 8*(-1)**n*n**2 - 21*(-1)**n*n + 18*(-1)**n]]"
+      ]
+     },
+     "execution_count": 58,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "table"
    ]
   },
   {