Ratio seems to determine relative error in each iteration

Author: hirish99

test/gridfree_taylor_recurrence.ipynb

@@ -9,7 +9,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -35,19 +35,18 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
     "var = _make_sympy_vec(\"x\", 2)\n",
     "s = sp.Function(\"s\")\n",
-    "n = sp.symbols(\"n\")\n",
-    "i = sp.symbols(\"i\")"
+    "n = sp.symbols(\"n\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -59,70 +58,81 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "var = _make_sympy_vec(\"x\", 2)\n",
+    "var_t = _make_sympy_vec(\"t\", 2)\n",
+    "g_x_y = sp.log(sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2))\n",
+    "derivs_laplace = [g_x_y.subs(var_t[0], 0).subs(var_t[1], 0).diff(var[0], i)\n",
+    "                    for i in range(8)]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/latex": [
-       "$\\displaystyle \\frac{2 s{\\left(1 \\right)}}{x_{1}^{2}}$"
+       "$\\displaystyle - \\frac{5040}{x_{1}^{8}}$"
       ],
       "text/plain": [
-       "2*s(1)/x1**2"
+       "-5040/x1**8"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 23,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "get_reindexed_and_center_origin_off_axis_recurrence(laplace2d)[2].subs(n, 3)"
+    "i=4\n",
+    "j=4\n",
+    "derivs_laplace[i].diff(var[0], j).subs(var[0], 0)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "(x0**3*(((-6*n + (n + 2)**2 - 4)*s(n - 3)/x1**2 + (-5*n + (n + 2)**2 - 4)*s(n - 1)/x1**2)*(-5*n + (n + 4)**2 - 14)/(6*x1**2) + (-6*n + (n + 4)**2 - 16)*s(n - 1)/(6*x1**2)) + x0**2*((-6*n + (n + 3)**2 - 10)*s(n - 2)/(2*x1**2) + (-5*n + (n + 3)**2 - 9)*s(n)/(2*x1**2)) + x0*((-6*n + (n + 2)**2 - 4)*s(n - 3)/x1**2 + (-5*n + (n + 2)**2 - 4)*s(n - 1)/x1**2) + (-6*n + (n + 1)**2 + 2)*s(n - 4)/x1**2 + (-5*n + (n + 1)**2 + 1)*s(n - 2)/x1**2,\n",
-       " 4)"
-      ]
-     },
-     "execution_count": 9,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "exp = get_off_axis_expression(helmholtz2d)\n",
-    "exp"
+    "k = 1\n",
+    "abs_dist = sp.sqrt((var[0]-var_t[0])**2 +\n",
+    "                       (var[1]-var_t[1])**2)\n",
+    "g_x_y = (1j/4) * hankel1(0, k * abs_dist)\n",
+    "derivs_helmholtz = [sp.diff(g_x_y,\n",
+    "                    var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)\n",
+    "                                            for i in range(6)]"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
-     "data": {
-      "text/latex": [
-       "$\\displaystyle \\frac{x_{0}^{3} \\left(- x_{1}^{2} \\left(6 n - \\left(n + 4\\right)^{2} + 16\\right) s{\\left(n - 1 \\right)} + \\left(\\left(5 n - \\left(n + 2\\right)^{2} + 4\\right) s{\\left(n - 1 \\right)} + \\left(6 n - \\left(n + 2\\right)^{2} + 4\\right) s{\\left(n - 3 \\right)}\\right) \\left(5 n - \\left(n + 4\\right)^{2} + 14\\right)\\right) + 3 x_{1}^{2} \\left(- x_{0}^{2} \\left(\\left(5 n - \\left(n + 3\\right)^{2} + 9\\right) s{\\left(n \\right)} + \\left(6 n - \\left(n + 3\\right)^{2} + 10\\right) s{\\left(n - 2 \\right)}\\right) - 2 x_{0} \\left(\\left(5 n - \\left(n + 2\\right)^{2} + 4\\right) s{\\left(n - 1 \\right)} + \\left(6 n - \\left(n + 2\\right)^{2} + 4\\right) s{\\left(n - 3 \\right)}\\right) + 2 \\left(- 6 n + \\left(n + 1\\right)^{2} + 2\\right) s{\\left(n - 4 \\right)} + 2 \\left(- 5 n + \\left(n + 1\\right)^{2} + 1\\right) s{\\left(n - 2 \\right)}\\right)}{6 x_{1}^{4}}$"
-      ],
-      "text/plain": [
-       "(x0**3*(-x1**2*(6*n - (n + 4)**2 + 16)*s(n - 1) + ((5*n - (n + 2)**2 + 4)*s(n - 1) + (6*n - (n + 2)**2 + 4)*s(n - 3))*(5*n - (n + 4)**2 + 14)) + 3*x1**2*(-x0**2*((5*n - (n + 3)**2 + 9)*s(n) + (6*n - (n + 3)**2 + 10)*s(n - 2)) - 2*x0*((5*n - (n + 2)**2 + 4)*s(n - 1) + (6*n - (n + 2)**2 + 4)*s(n - 3)) + 2*(-6*n + (n + 1)**2 + 2)*s(n - 4) + 2*(-5*n + (n + 1)**2 + 1)*s(n - 2)))/(6*x1**4)"
-      ]
-     },
-     "execution_count": 11,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "NameError",
+     "evalue": "name 'derivs_helmholtz' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mderivs_helmholtz\u001b[49m[\u001b[38;5;241m4\u001b[39m]\u001b[38;5;241m.\u001b[39mdiff(var[\u001b[38;5;241m1\u001b[39m], \u001b[38;5;241m0\u001b[39m)\u001b[38;5;241m.\u001b[39msubs(var[\u001b[38;5;241m1\u001b[39m], \u001b[38;5;241m0\u001b[39m)\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'derivs_helmholtz' is not defined"
+     ]
     }
    ],
    "source": [
-    "exp[0].simplify()"
+    "derivs_helmholtz[4].diff(var[1], 0).subs(var[1], 0)"
    ]
   },
   {

test/investigate_normal.py

@@ -0,0 +1,148 @@
+# %%
+from sumpy.recurrence import _make_sympy_vec, get_reindexed_and_center_origin_on_axis_recurrence
+
+from sumpy.expansion.diff_op import (
+    laplacian,
+    make_identity_diff_op,
+)
+
+
+import sympy as sp
+from sympy import hankel1
+
+import numpy as np
+
+import math
+
+import matplotlib.pyplot as plt
+from matplotlib import cm, ticker
+
+# %%
+w = make_identity_diff_op(2)
+laplace2d = laplacian(w)
+n_init_lap, order_lap, recur_laplace = get_reindexed_and_center_origin_on_axis_recurrence(laplace2d)
+
+w = make_identity_diff_op(2)
+helmholtz2d = laplacian(w) + w
+n_init_helm, order_helm, recur_helmholtz = get_reindexed_and_center_origin_on_axis_recurrence(helmholtz2d)
+
+# %%
+var = _make_sympy_vec("x", 2)
+rct = sp.symbols("r_{ct}")
+s = sp.Function("s")
+n = sp.symbols("n")
+
+# %%
+def compute_derivatives(p):
+    var = _make_sympy_vec("x", 2)
+    var_t = _make_sympy_vec("t", 2)
+    g_x_y = sp.log(sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2))
+    derivs = [sp.diff(g_x_y,
+                        var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)
+                        for i in range(p)]
+    return derivs
+l_max = 10
+derivs_laplace = compute_derivatives(l_max)
+
+# %%
+def compute_derivatives_h2d(p):
+    k = 1
+    var = _make_sympy_vec("x", 2)
+    var_t = _make_sympy_vec("t", 2)
+    abs_dist = sp.sqrt((var[0]-var_t[0])**2 +
+                        (var[1]-var_t[1])**2)
+    g_x_y = (1j/4) * hankel1(0, k * abs_dist)
+    derivs_helmholtz = [sp.diff(g_x_y,
+                        var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0)
+                                                for i in range(p)]
+    return derivs_helmholtz
+h_max = 8
+#derivs_helmholtz = compute_derivatives_h2d(h_max)
+
+# %%
+def evaluate_recurrence_lamb(coord_dict, recur, p, derivs_list, n_initial, n_order):
+    s = sp.Function("s")
+    subs_dict = {}
+    for i in range(n_initial-n_order, 0):
+        subs_dict[s(i)] = 0
+    for i in range(n_initial):
+        subs_dict[s(i)] = derivs_list[i].subs(coord_dict)
+    var = _make_sympy_vec("x", 2)
+    for i in range(n_initial, p):
+        exp = recur.subs(n, i)
+        f = sp.lambdify([var[0], var[1]] + [s(i-(1+k)) for k in range(n_order-1)], exp)
+        subs_dict[s(i)] = f(*([coord_dict[var[0]], coord_dict[var[1]]] + [subs_dict[s(i-(1+k))] for k in range(n_order-1)]))
+    for i in range(n_initial-n_order, 0):
+        subs_dict.pop(s(i))
+    return np.array(list(subs_dict.values()))
+
+# %%
+def evaluate_true(coord_dict, p, derivs_list):
+    retMe = []
+    for i in range(p):
+        exp = derivs_list[i]
+        f = sp.lambdify(var, exp)
+        retMe.append(f(coord_dict[var[0]], coord_dict[var[1]]))
+    return np.array(retMe)
+
+# %%
+def compute_error_coord(recur, loc, order, derivs_list, n_initial, n_order):
+    var = _make_sympy_vec("x", 2)
+    coord_dict = {var[0]: loc[0], var[1]: loc[1]}
+
+    exp = evaluate_recurrence_lamb(coord_dict, recur, order+1, derivs_list, n_initial, n_order)[order].evalf()
+    
+    true = derivs_list[order].subs(coord_dict).evalf()
+
+    return (np.abs(exp-true)/np.abs(true))
+
+# %%
+def generate_error_grid(res, order_plot, recur, derivs, n_initial, n_order):
+    x_grid = [10**(pw) for pw in np.linspace(-8, 0, res)]
+    y_grid = [10**(pw) for pw in np.linspace(-8, 0, res)]
+    res=len(x_grid)
+    plot_me = np.empty((res, res))
+    for i in range(res):
+        for j in range(res):
+            if abs(x_grid[i]) == abs(y_grid[j]):
+                plot_me[i, j] = 1e-16
+            else:
+                plot_me[i,j] = compute_error_coord(recur, np.array([x_grid[i],y_grid[j]]), order_plot, derivs, n_initial, n_order)
+                if plot_me[i,j] == 0:
+                    plot_me[i, j] = 1e-16
+    return x_grid, y_grid, plot_me
+
+# %%
+order_plot = 5
+#x_grid, y_grid, plot_me_hem = generate_error_grid(res=5, order_plot=order_plot, recur=recur_helmholtz, derivs=derivs_helmholtz, n_initial=n_init_helm, n_order=order_helm)
+x_grid, y_grid, plot_me_lap  = generate_error_grid(res=10, order_plot=order_plot, recur=recur_laplace, derivs=derivs_laplace, n_initial=n_init_lap, n_order=order_lap)
+plot_me_hem = plot_me_lap
+      
+                
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 8))
+cs = ax1.contourf(x_grid, y_grid, plot_me_hem.T, locator=ticker.LogLocator(), cmap=cm.PuBu_r)
+cbar = fig.colorbar(cs)
+
+cs = ax2.contourf(x_grid, y_grid, plot_me_lap.T, locator=ticker.LogLocator(), cmap=cm.PuBu_r)
+cbar = fig.colorbar(cs)
+ax1.set_xscale('log')
+ax1.set_yscale('log')
+ax1.set_xlabel("source x-coord", fontsize=15)
+ax1.set_ylabel("source y-coord", fontsize=15)
+
+
+ax2.set_xscale('log')
+ax2.set_yscale('log')
+ax2.set_xlabel("source x-coord", fontsize=15)
+ax2.set_ylabel("source y-coord", fontsize=15)
+
+ax1.set_title("Helmholtz recurrence relative error for order = "+str(order_plot), fontsize=15)
+ax2.set_title("Laplace recurrence relative error for order = "+str(order_plot), fontsize=15)
+
+fig.savefig('order'+str(order_plot))
+plt.show()
+
+# %%
+
+
+