First_Order_SBP_Operators.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f697c086",
   "metadata": {},
   "source": [
    "# SBP Operators\n",
    "#### References:\n",
    "\n",
    "Carpenter, M. H., Parsani, M., Fisher, T. C., & Nielsen, E. J. (2015). Entropy Stable Staggered Grid Spectral Collocation for the Burgers' and Compressible Navier-Stokes Equations (No. NF1676L-22150).\n",
    "\n",
    "Canuto, C., Hussaini, M. Y., Quarteroni, A., & Zang, T. A. (2007). Spectral methods: fundamentals in single domains. Springer Science & Business Media.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c497864",
   "metadata": {},
   "source": [
    "## Computing Lagrange basis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eee9b848",
   "metadata": {},
   "source": [
    "We import some useful modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f4555e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26c76edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from basis_nodes import generate_lagrange_poly, generate_LGL_points"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcbcedc9",
   "metadata": {},
   "source": [
    "The function generate_lagrange_poly(j, x_nodes) will return a Lagrange polynomial and its derivative:\n",
    "\n",
    "$$\n",
    "L_j(x)=\\prod_{k \\neq j} \\frac{x-x_k}{x_j-x_k},\n",
    "$$\n",
    "and\n",
    "$$\n",
    "L_j^{\\prime}(x)=\\frac{\\sum_{k=0, k \\neq j}^n \\prod_{l=0, l \\neq k, l \\neq j}^n\\left(x-x_l\\right)}{\\prod_{k=0, k \\neq j}^n\\left(x_j-x_k\\right)}.\n",
    "$$\n",
    "\n",
    "A set of nodes should be provided, it is desired that these be the Legendre-Gauss-Lobatto (LGL) nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15aa99ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = np.array([0,1,2,3,4,5,6])\n",
    "idx = 4\n",
    "poly, polyp = generate_lagrange_poly(idx,nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aa355f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "vpoly = np.vectorize(poly)\n",
    "vpolyp = np.vectorize(polyp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cea8e8fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,6,50)\n",
    "\n",
    "fig = plt.figure(dpi=100, figsize=(8,4))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_ylim([-5,5])\n",
    "ax.plot(x,np.zeros(len(x)),color='k', lw=1)\n",
    "ax.plot(x, vpoly(x), label=f'$L_{idx}(x)$')\n",
    "ax.plot(x, vpolyp(x),color='g', label=f'$L_{idx}\\'(x)$')\n",
    "ax.scatter(nodes,vpoly(nodes),marker='x', color='r')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd7b17e2",
   "metadata": {},
   "source": [
    "## Computing LG and LGL points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d71faed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "degree_basis = 10 #degree_basis = n_points-1\n",
    "(p_Legn, p_Legn_prime,\n",
    "p_Legn1, p_Legn1_prime,\n",
    "x_LG_unsorted, w_LG_unsorted,\n",
    "x_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
    "\n",
    "#####\n",
    "#The root finder returns unsorted roots,\n",
    "#We fix that in these lines\n",
    "#####\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LG = np.argsort(x_LG_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "x_LG = x_LG_unsorted[sort_idxs_LG]\n",
    "w_LG = w_LG_unsorted[sort_idxs_LG]\n",
    "\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LGL = np.argsort(x_LGL_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "x_LGL = x_LGL_unsorted[sort_idxs_LGL]\n",
    "w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
    "\n",
    "x = np.linspace(-1,1,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fff31b9",
   "metadata": {},
   "source": [
    "The Legendre polynomials are computed as eigenfunctions of the Sturm-Liouville problem (Canuto et al., 2006):\n",
    "\n",
    "$$\n",
    "((1-x^2) L_k'(x))'+k(k+1) L_k(x) = 0.\n",
    "$$\n",
    "\n",
    "The Legendre-Gauss nodes and weights are computed as\n",
    "$$\n",
    "\\begin{gathered}\n",
    "x_j(j=0, \\ldots, N) \\text { zeros of } L_{N+1} ; \\\\\n",
    "w_j=\\frac{2}{\\left(1-x_j^2\\right)\\left[L_{N+1}^{\\prime}\\left(x_j\\right)\\right]^2}, \\quad j=0, \\ldots, N\n",
    "\\end{gathered}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8bf54bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(dpi=100, figsize=(8,3))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(x,np.zeros(len(x)),color='k', lw=1)\n",
    "ax.plot(x, p_Legn1(x), label=f'$L(x)$')\n",
    "ax.scatter(x_LG,np.zeros(len(x_LG)),marker='x', color='r', label='LG points')\n",
    "#ax.scatter(x_LG,w_LG,marker='x', color='g', label='LGL weights')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f20a5d1",
   "metadata": {},
   "source": [
    "The Legendre-Gauss-Lobatto nodes and weights are computed as\n",
    "\n",
    "$$\n",
    "\\begin{gathered}\n",
    "x_0=-1, x_N=1, x_j(j=1, \\ldots, N-1) \\text { zeros of } L_N^{\\prime} ; \\\\\n",
    "w_j=\\frac{2}{N(N+1)} \\frac{1}{\\left[L_N\\left(x_j\\right)\\right]^2}, \\quad j=0, \\ldots, N .\n",
    "\\end{gathered}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4febe4c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#tol=1.e-14\n",
    "#findiff=(p_Legn(x+tol)-p_Legn(x-tol))/tol\n",
    "fig = plt.figure(dpi=100, figsize=(8,3))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.plot(x,np.zeros(len(x)),color='k', lw=1)\n",
    "ax.plot(x, p_Legn_prime(x), label=f'$L\\'(x)$')\n",
    "#ax.plot(x,findiff, color='g')\n",
    "ax.scatter(x_LGL,np.zeros(len(x_LGL)),marker='x', color='r', label='LGL points')\n",
    "#ax.scatter(x_LGL,w_LGL,marker='x', color='g', label='LGL weights')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9246ea0",
   "metadata": {},
   "source": [
    "Validating nodes and weights with some simple integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40795bd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "g = lambda x: (1./3.)*x**3\n",
    "num_int_LG=np.sum(f(x_LG)*w_LG)\n",
    "num_int_LGL=np.sum(f(x_LGL)*w_LGL)\n",
    "exact_int = g(1)-g(-1)\n",
    "print(f\"LG integration: {num_int_LG}\")\n",
    "print(f\"LGL integration:{num_int_LGL}\")\n",
    "print(f\"Exact integral: {exact_int}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e94c299",
   "metadata": {},
   "source": [
    "## Computing first order SBP matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd474626",
   "metadata": {},
   "source": [
    "### First Order Differential Operator\n",
    "$$\n",
    "\\mathcal{D}=\\left(d_{i j}\\right)=\\left(\\frac{d L_j}{d x}\\left(x_i\\right)\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2db3c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from SBP_matrices import first_order_D, first_order_P_Q\n",
    "\n",
    "degree_basis = 2 #degree_basis = n_nodes-1\n",
    "(p_Legn, p_Legn_prime,\n",
    "p_Legn1, p_Legn1_prime,\n",
    "x_LG_unsorted, w_LG_unsorted,\n",
    "x_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
    "\n",
    "#####\n",
    "#The root finder returns unsorted roots,\n",
    "#We fix that in these lines\n",
    "#####\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LG = np.argsort(x_LG_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "x_LG = x_LG_unsorted[sort_idxs_LG]\n",
    "w_LG = w_LG_unsorted[sort_idxs_LG]\n",
    "\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LGL = np.argsort(x_LGL_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "x_LGL = x_LGL_unsorted[sort_idxs_LGL]\n",
    "w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
    "\n",
    "n_nodes = len(x_LGL)\n",
    "x = np.linspace(-0.5,0.5,n_nodes)\n",
    "D_LG = first_order_D(x_nodes=x_LG)\n",
    "D_LGL = first_order_D(x_nodes=x_LGL)\n",
    "D = first_order_D(x_nodes=x)\n",
    "print(\"D_LG:\")\n",
    "print(np.round(D_LG,2))\n",
    "print(\"D_LGL:\")\n",
    "print(np.round(D_LGL,2))\n",
    "print(\"D:\")\n",
    "print(np.round(D,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c42c58a3",
   "metadata": {},
   "source": [
    "### Computing P and Q\n",
    "\n",
    "$$\n",
    "\\mathcal{P}=\\sum_{\\ell} \\mathbf{L}\\left(\\eta_l ; \\mathbf{x}\\right) \\mathbf{L}\\left(\\eta_l ; \\mathbf{x}\\right)^{\\top} \\omega_{\\ell} \\quad, \\quad \\mathcal{Q}=\\sum_{\\ell} \\mathbf{L}\\left(\\eta_l ; \\mathbf{x}\\right) \\frac{d \\mathbf{L}}{d x}\\left(\\eta_l ; \\mathbf{x}\\right)^{\\top} \\omega_{\\ell}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3f14561",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_LG, Q_LG = first_order_P_Q(x_Lagrange_nodes=x_LG, x_abcissae=x_LG, w_abcissae=w_LG)\n",
    "P_LGL, Q_LGL = first_order_P_Q(x_Lagrange_nodes=x_LGL, x_abcissae=x_LGL, w_abcissae=w_LGL)\n",
    "P, Q = first_order_P_Q(x_Lagrange_nodes=x, x_abcissae=x_LGL, w_abcissae=w_LGL)\n",
    "print(\"P_LG:\")\n",
    "print(np.round(P_LG,2))\n",
    "print(\"P_LGL:\")\n",
    "print(np.round(P_LGL,2))\n",
    "print(\"P:\")\n",
    "print(np.round(P,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80f8fbcf",
   "metadata": {},
   "source": [
    "### Validating our computations using Q\n",
    "We compute $Q$ through $Q=PQ$, and compare it with the matrices we obtained before\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fae9e851",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"P_LG @ D_LG:\")\n",
    "print(np.round(P_LG@D_LG,2))\n",
    "print(\"P_LGL @ D_LG:\")\n",
    "print(np.round(P_LGL@D_LGL,2))\n",
    "print(\"P @ D:\")\n",
    "print(np.round(P@D,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f332f16",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Q_LG:\")\n",
    "print(np.round(Q_LG,2))\n",
    "print(\"Q_LGL:\")\n",
    "print(np.round(Q_LGL,2))\n",
    "print(\"Q:\")\n",
    "print(np.round(Q,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99a460ff",
   "metadata": {},
   "source": [
    "### Verifying SBP properties\n",
    "From inspection we can see that $P$ is symmetric for the different Lagrange bases, lets check if it's PD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "052e7059",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_pos_def(A):\n",
    "    if np.array_equal(A, A.T):\n",
    "        try:\n",
    "            np.linalg.cholesky(A)\n",
    "            return True\n",
    "        except np.linalg.LinAlgError:\n",
    "            return False\n",
    "    else:\n",
    "        return False\n",
    "print(f\"P_LG SPD?: {is_pos_def(P_LG)}\")\n",
    "print(f\"P_LGL SPD?: {is_pos_def(P_LGL)}\")\n",
    "print(f\"P SPD?: {is_pos_def(P)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bda43954",
   "metadata": {},
   "source": [
    "Now we check if\n",
    "$$\n",
    "    Q^T = B- Q. \\qquad B=\\text{Diag}(-1,0,...,0,1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7426fb60",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "B=np.zeros((len(P),len(P))); B[0,0]=-1; B[-1,-1]=1\n",
    "print(f\"Q_LG^T - (B-Q_LG):\\n\", np.round(np.transpose(Q_LG)-(B-Q_LG),2))\n",
    "print(f\"Q_LGL^T - (B-Q_LGL):\\n\", np.round(np.transpose(Q_LGL)-(B-Q_LGL),2))\n",
    "print(f\"Q^T - (B-Q):\\n\", np.round(np.transpose(Q)-(B-Q),2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14551135",
   "metadata": {},
   "source": [
    "Apart from having a diagonal matrix $P_{LGL}$, the Lagrange basis produced by the LGL nodes gives us matrices $P_{LGL}$ and $Q_{LGL}$ that satisfy the SBP convention."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addba995",
   "metadata": {
    "hideCode": true
   },
   "source": [
    "### Testing action of $D$ and  $D^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bcf3b51",
   "metadata": {},
   "outputs": [],
   "source": [
    "degree_basis = 4 #degree_basis = n_nodes-1\n",
    "(p_Legn, p_Legn_prime,\n",
    "p_Legn1, p_Legn1_prime,\n",
    "x_LG_unsorted, w_LG_unsorted,\n",
    "x_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LGL = np.argsort(x_LGL_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "x_LGL = x_LGL_unsorted[sort_idxs_LGL]\n",
    "w_LGL = w_LGL_unsorted[sort_idxs_LGL]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27ee2d8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "D_LGL = first_order_D(x_nodes=x_LGL)\n",
    "P_LGL, Q_LGL = first_order_P_Q(x_Lagrange_nodes=x_LGL, x_abcissae=x_LGL, w_abcissae=w_LGL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c019acb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ftest = lambda x: 2*x**2+x**3\n",
    "plt.plot(x_LGL,4*x_LGL+3*x_LGL**2)\n",
    "plt.plot(x_LGL,D_LGL@ftest(x_LGL))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b96145a",
   "metadata": {},
   "source": [
    "## LGL nodes and weights, and SBP operators mapped to a new interval\n",
    "\n",
    "# Warning: Computing Lagrange basis at each physical element leads to instablities, it's better to work on the reference element and just keep track of the mappings from computational to physical domain."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba077555",
   "metadata": {},
   "source": [
    "It is often desired to work in a different interval than $[-1,1]$, say $[x_\\min, x_\\max]$. \n",
    "We can map the LGL nodes using the affine mapping:\n",
    "\n",
    "$$\n",
    "x = a \\xi + b,\\quad \\xi\\in[-1,1],\\; x\\in [x_\\min,x_\\max],\n",
    "$$\n",
    "\n",
    "where\n",
    "$$\n",
    "a = \\frac{x_\\max-x_\\min}{2}, \\quad b=\\frac{x_\\min+x_\\max}{2}.\n",
    "$$\n",
    "\n",
    "Furthermore, we can approximate integrals as\n",
    "\n",
    "$$\n",
    "\\int_{x_\\min}^{x_\\max} g(x)dx = a\\int_{-1}^1 g(a\\xi+b)d\\xi \\approx a\\sum_{i=1}^N w_i g(a\\xi_i+b) =\\sum_{i=1}^N a w_i g(x_i).\n",
    "$$\n",
    "\n",
    "Therefore, to obtain weights for the new mapped abcissae, it suffices to multiply the LGL weights by $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b1296b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "degree_basis = 4 #degree_basis = n_nodes-1\n",
    "(_, _,\n",
    "_, _,\n",
    "_, _,\n",
    "xi_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
    "\n",
    "#######################\n",
    "#######################\n",
    "#The root finder returns unsorted roots,\n",
    "#We fix that in these lines\n",
    "\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LGL = np.argsort(xi_LGL_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "xi_LGL = xi_LGL_unsorted[sort_idxs_LGL]\n",
    "w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
    "\n",
    "xmin = -3; xmax=3\n",
    "a = (xmax-xmin)/2.; b = (xmax+xmin)/2.\n",
    "x_mapped = a*xi_LGL+b\n",
    "w_mapped = a*w_LGL\n",
    "#######################\n",
    "#######################\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89b8b473",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"x={x_mapped}\")\n",
    "print(f\"Sum of weights\")\n",
    "print(np.sum(w_LGL))\n",
    "print(np.sum(w_mapped))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b660f3",
   "metadata": {},
   "source": [
    "Since the Lagrange basis and, therefore, the first order differential operator just require a set of nodes for computation, they are computed exactly as before.\n",
    "However, special care should be taken for the computation of $P$ and $Q$.\n",
    "### Lagrange basis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9eb7843",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_support = 3 #Which node is going to be 1\n",
    "poly, polyp = generate_lagrange_poly(idx_support,x_mapped)\n",
    "x = np.linspace(-3.5,3.5,50)\n",
    "vpoly = np.vectorize(poly)\n",
    "vpolyp = np.vectorize(polyp)\n",
    "fig = plt.figure(dpi=100, figsize=(8,4))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_ylim([-5,5])\n",
    "ax.plot(x,np.zeros(len(x)),color='k', lw=1)\n",
    "ax.plot(x, vpoly(x), label=f'$L_{idx_support}(x)$')\n",
    "ax.plot(x, vpolyp(x),color='g', label=f'$L_{idx_support}\\'(x)$')\n",
    "ax.scatter(x_mapped,vpoly(x_mapped),marker='x', color='r')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8f02bdd",
   "metadata": {},
   "source": [
    "### First Order Differential Operator D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4507dee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "D_mapped = first_order_D(x_nodes=x_mapped)\n",
    "print(np.round(D_mapped,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8b92631",
   "metadata": {},
   "source": [
    "### Computing P and Q\n",
    "\n",
    "The computation of P and Q require the evaluation of the Lagrange basis (obtained with the collocation points $\\mathbf{x}$) in the LGL nodes $\\eta$. Therefore, the mapped abcissae and weights should be used to generate and evaluate the Lagrange polynomials:\n",
    "$$\n",
    "\\mathcal{P}= a\\sum_{\\ell} \\mathbf{L}\\left(x(\\eta_l) ; \\mathbf{x}\\right) \\mathbf{L}\\left(x(\\eta_l) ; \\mathbf{x}\\right)^{\\top} \\omega_{\\ell} \\quad, \\quad \n",
    "\\mathcal{Q}=a\\sum_{\\ell} \\mathbf{L}\\left(x(\\eta_l) ; \\mathbf{x}\\right) \\frac{d \\mathbf{L}}{d x}\\left(x(\\eta_l) ; \\mathbf{x}\\right)^{\\top} \\omega_{\\ell}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f690af2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "P, Q = first_order_P_Q(x_Lagrange_nodes=x_mapped, x_abcissae=x_mapped, w_abcissae=w_mapped)\n",
    "print(\"P:\")\n",
    "print(np.round(P,2))\n",
    "print(\"Q:\")\n",
    "print(np.round(Q,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "585fcfa2",
   "metadata": {},
   "source": [
    "### Verifying SBP properties\n",
    "We check if P is SPD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3b4dc89",
   "metadata": {},
   "outputs": [],
   "source": [
    "is_pos_def(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ccd0d8b",
   "metadata": {},
   "source": [
    "Now we check if\n",
    "$$\n",
    "    Q^T = B- Q. \\qquad B=\\text{Diag}(-1,0,...,0,1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46f1b6a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "B=np.zeros((len(P),len(P))); B[0,0]=-1; B[-1,-1]=1\n",
    "print(f\"Q^T - (B-Q):\\n\", np.round(np.transpose(Q)-(B-Q),2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.8.9",
   "language": "python",
   "name": "pyenv_3.8.9"
  },
  "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.8.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}