{ "cells": [ { "cell_type": "markdown", "id": "2bb2025a", "metadata": { "editable": true }, "source": [ "\n", "\n", "# Demo - Sparse Chebyshev-Petrov-Galerkin methods for differentiation\n", "**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n", "\n", "Date: **October 26, 2021**\n", "\n", "**Summary.** This demo explores how to use sparse Chebyshev-Petrov-Galerkin methods for finding Chebyshev coefficients of\n", "the derivatives of smooth functions. We will compare the methods to the more commonly adopted\n", "recursion methods that are found in most spectral textbooks." ] }, { "cell_type": "markdown", "id": "4838d7e1", "metadata": { "editable": true }, "source": [ "## Introduction\n", "\n", "The Chebyshev polynomials of the first kind can be defined as" ] }, { "cell_type": "markdown", "id": "1af52674", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eq:chebTU} \\tag{1}\n", " T_k(x) = \\cos(k\\theta),\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "747a43e9", "metadata": { "editable": true }, "source": [ "where $\\theta = \\cos^{-1} x$, $k$ is a positive integer and $x \\in [-1, 1]$.\n", "The Chebyshev polynomials span the discrete space $S_N = \\text{span}\\{T_k\\}_{k=0}^{N-1}$,\n", "and a function $u(x)$ can be approximated in this space as" ] }, { "cell_type": "markdown", "id": "65fce623", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u_N(x) = \\sum_{k=0}^{N-1} \\hat{u}_k T_k(x).\n", "\\label{eq:uT} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1f1b1a55", "metadata": { "editable": true }, "source": [ "Consider the expansion of the function $u(x)=\\sin(\\pi x)$, created in `shenfun` as" ] }, { "cell_type": "code", "execution_count": 1, "id": "e8d08068", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:56.145114Z", "iopub.status.busy": "2024-02-19T13:48:56.144839Z", "iopub.status.idle": "2024-02-19T13:48:57.016848Z", "shell.execute_reply": "2024-02-19T13:48:57.016336Z" } }, "outputs": [ { "data": { "text/plain": [ "Function([-4.85722573e-17, 5.69230686e-01, -1.38366483e-16,\n", " -6.66916672e-01, 7.03290092e-17, 1.04282369e-01,\n", " -5.17536910e-17, -6.84063354e-03, -2.94392336e-17,\n", " 2.50006885e-04, 8.57780639e-18, -5.85024831e-06,\n", " -6.09960263e-17, 9.53478051e-08, -1.04132138e-16,\n", " -1.15621280e-09])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from shenfun import *\n", "import sympy as sp\n", "x = sp.Symbol('x')\n", "ue = sp.sin(sp.pi*x)\n", "N = 16\n", "SN = FunctionSpace(N, 'C')\n", "uN = Function(SN, buffer=ue)\n", "uN" ] }, { "cell_type": "markdown", "id": "bb86e5fd", "metadata": { "editable": true }, "source": [ "The Python Function `uN` represents the expansion ([2](#eq:uT)), and the printed\n", "values represent $\\boldsymbol{\\hat{u}} = \\{\\hat{u}_k\\}_{k=0}^{N-1}$. The expansion is fairly well resolved since\n", "the highest values of $\\{\\hat{u}_k\\}_{k=0}^{N-1}$ approach 0.\n", "Note that the coefficients obtained are based on interpolation at\n", "quadrature points and they do not agree completely with the coefficients truncated from an\n", "infinite series $u(x) = \\sum_{k=0}^{\\infty} \\hat{u}_k T_k$. The obtained series is\n", "often denoted as $u_N(x) = I_N u(x)$, where $I_N$ is an interpolation operator.\n", "Under the hood the coefficients are found by projection using quadrature for the integrals:\n", "find $u_N \\in S_N$ such that" ] }, { "cell_type": "markdown", "id": "563fe691", "metadata": { "editable": true }, "source": [ "$$\n", "(u_N-u, v)_{\\omega^{-1/2}} = 0, \\quad \\forall v \\in S_N,\n", "$$" ] }, { "cell_type": "markdown", "id": "846b9e53", "metadata": { "editable": true }, "source": [ "where $\\omega = (1-x^2)$ and the scalar product notation\n", "$(a, b)_{\\omega^{-1/2}} = \\sum_{j=0}^{N-1} a(x_j)b(x_j)\\omega_j \\approx \\int_{-1}^{1} a(x)b(x) \\omega(x)^{-1/2} dx$,\n", "where $\\{\\omega_j\\}_{j=0}^{N-1}$ are the quadrature weights. The quadrature approach ensures\n", "that $u(x_j) = u_N(x_j)$ for all quadrature points $\\{x_j\\}_{j=0}^{N-1}$.\n", "In `shenfun` we compute the following under the hood: insert for $u_N = \\sum_{j=0}^{N-1} \\hat{u}_j T_j$,\n", "$u=\\sin(\\pi x)$ and $v = T_k$ to get" ] }, { "cell_type": "markdown", "id": "833df5cb", "metadata": { "editable": true }, "source": [ "$$\n", "\\sum_{j=0}^{N-1}(T_j, T_k)_{\\omega^{-1/2}} \\hat{u}_j = (\\sin(\\pi x), T_k)_{\\omega^{-1/2}},\n", "$$" ] }, { "cell_type": "markdown", "id": "54d5f8e3", "metadata": { "editable": true }, "source": [ "This has now become a linear algebra problem, and we recognise the matrix $d^{(0)}_{kj} = (T_j, T_k)_{\\omega^{-1/2}}=c_k \\pi /2 \\delta_{kj}$,\n", "where $\\delta_{kj}$ is the Kronecker delta function, and $c_0=2$ and $c_k=1$ for $k>0$.\n", "The problem is solved trivially since $d^{(0)}_{kj}$ is diagonal,\n", "and thus" ] }, { "cell_type": "markdown", "id": "1ad4aac0", "metadata": { "editable": true }, "source": [ "$$\n", "\\hat{u}_k = \\frac{2}{c_k \\pi} (\\sin(\\pi x), T_k)_{\\omega^{-1/2}}, \\quad \\forall \\, k\\in I^N,\n", "$$" ] }, { "cell_type": "markdown", "id": "ced04b4b", "metadata": { "editable": true }, "source": [ "where $I^N = \\{0, 1, \\ldots, N-1\\}$.\n", "We can compare this to the truncated coefficients, where the integral $(\\sin(\\pi x), T_k)_{\\omega^{-1/2}}$\n", "is computed with high precision. To this end we could use adaptive quadrature, or symbolic integration\n", "with sympy, but it is sufficient to use a large enough number of polynomials to fully resolve the\n", "function. Below we find this number to be 22 and we see that the absolute error in\n", "the highest wavenumber $\\hat{u}_{15} \\approx 10^{-11}$." ] }, { "cell_type": "code", "execution_count": 2, "id": "cd4b2ee5", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.020267Z", "iopub.status.busy": "2024-02-19T13:48:57.019937Z", "iopub.status.idle": "2024-02-19T13:48:57.032588Z", "shell.execute_reply": "2024-02-19T13:48:57.032143Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 8.20051098e-18 1.11022302e-16 1.22890593e-17 -1.11022302e-16\n", " -4.29912860e-18 5.55111512e-17 8.53325551e-17 -2.44596010e-16\n", " 1.85259266e-17 -6.84673672e-17 6.31265804e-17 4.99434343e-16\n", " 1.63192838e-16 -7.76732293e-14 1.21113576e-16 1.05743425e-11]\n", "22\n" ] } ], "source": [ "SM = FunctionSpace(0, 'C')\n", "uM = Function(SM, buffer=ue, abstol=1e-16, reltol=1e-16)\n", "print(uM[:N] - uN[:N])\n", "print(len(uM))" ] }, { "cell_type": "markdown", "id": "bef631d0", "metadata": { "editable": true }, "source": [ "## Differentiation\n", "\n", "Let us now consider the $n$'th derivative of $u(x)$ instead, denoted here as $u^{(n)}$.\n", "We want to find $u^{(n)}$ in the space $S_N$, which means that we want to obtain\n", "the following expansion" ] }, { "cell_type": "markdown", "id": "2c78531f", "metadata": { "editable": true }, "source": [ "$$\n", "u_N^{(n)} = \\sum_{k=0}^{N-1} \\hat{u}^{(n)}_k T_k.\n", "$$" ] }, { "cell_type": "markdown", "id": "6bacc81d", "metadata": { "editable": true }, "source": [ "First note that this expansion is not the same as the derivative of\n", "the previously found $u_N$, which is" ] }, { "cell_type": "markdown", "id": "aa75f2fd", "metadata": { "editable": true }, "source": [ "$$\n", "(u_N)^{(n)} = \\sum_{k=0}^{N-1} \\hat{u}_k T^{(n)}_k,\n", "$$" ] }, { "cell_type": "markdown", "id": "311527fa", "metadata": { "editable": true }, "source": [ "where $T^{(n)}_k$ is the $n$'th derivative of $T_k$, a polynomial of order $k-n$.\n", "We again use projection to find $u_N^{(n)} \\in S_N$ such that" ] }, { "cell_type": "markdown", "id": "5bf2f5d1", "metadata": { "editable": true }, "source": [ "$$\n", "(u_N^{(n)}-u^{(n)}, v)_{\\omega^{-1/2}} = 0, \\quad \\forall v \\in S_N.\n", "$$" ] }, { "cell_type": "markdown", "id": "f7e7d354", "metadata": { "editable": true }, "source": [ "Inserting for $u_N^{(n)}$ and $u^{(n)} = (u_N)^{(n)}$ we get" ] }, { "cell_type": "markdown", "id": "45bc180f", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\sum_{j=0}^{N-1}(T_j, T_k)_{\\omega^{-1/2}} \\hat{u}_j^{(n)} = (T_j^{(n)}, T_k)_{\\omega^{-1/2}} \\hat{u}_j, \n", "\\label{_auto1} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "bd058277", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " \\sum_{j=0}^{N-1} d^{(0)}_{kj} \\hat{u}_j^{(n)} = \\sum_{j=0}^{N-1} d^{(n)}_{kj} \\hat{u}_j,\n", "\\label{_auto2} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1dc093f5", "metadata": { "editable": true }, "source": [ "where $d^{(n)}_{kj} = (T_j^{(n)}, T_k)_{\\omega^{-1/2}}$.\n", "We compute $\\hat{u}_k^{(n)}$ by inverting the diagonal $d^{(0)}_{kj}$" ] }, { "cell_type": "markdown", "id": "f6025179", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\hat{u}_k^{(n)} = \\frac{2}{c_k \\pi} \\sum_{j=0}^{N-1} d^{(n)}_{kj} \\hat{u}_j, \\quad \\forall \\, k \\in I^{N}.\n", "\\label{eq:fhat} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "7eea7af8", "metadata": { "editable": true }, "source": [ "The matrix $d^{(n)}_{kj}$ is upper triangular, and the last $n$ rows are zero. Since $d^{(n)}_{kj}$ is\n", "dense the matrix vector product $\\sum_{j=0}^{N-1} d^{(n)}_{kj} \\hat{u}_j$ is costly\n", "and also susceptible to roundoff errors if the structure of the matrix is\n", "not taken advantage of. But computing it in shenfun\n", "is straightforward, for $n=1$ and $2$:" ] }, { "cell_type": "code", "execution_count": 3, "id": "4f6c7409", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.035793Z", "iopub.status.busy": "2024-02-19T13:48:57.035605Z", "iopub.status.idle": "2024-02-19T13:48:57.044542Z", "shell.execute_reply": "2024-02-19T13:48:57.044117Z" } }, "outputs": [ { "data": { "text/plain": [ "Function([-9.55804991e-01, -5.29095427e-15, -3.05007135e+00,\n", " -4.73748833e-15, 9.51428681e-01, -5.30012041e-15,\n", " -9.13950067e-02, -4.67907612e-15, 4.37386282e-03,\n", " -4.20804838e-15, -1.26261106e-04, -4.37960451e-15,\n", " 2.44435655e-06, -2.91569988e-15, -3.46863840e-08,\n", " 0.00000000e+00])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "uN1 = project(Dx(uN, 0, 1), SN)\n", "uN2 = project(Dx(uN, 0, 2), SN)\n", "uN1" ] }, { "cell_type": "markdown", "id": "ed23d795", "metadata": { "editable": true }, "source": [ "where `uN1` $=u_N^{(1)}$ and `uN2` $=u_N^{(2)}$.\n", "Alternatively, doing all the work that goes on under the hood" ] }, { "cell_type": "code", "execution_count": 4, "id": "85617695", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.046970Z", "iopub.status.busy": "2024-02-19T13:48:57.046775Z", "iopub.status.idle": "2024-02-19T13:48:57.053802Z", "shell.execute_reply": "2024-02-19T13:48:57.053104Z" } }, "outputs": [ { "data": { "text/plain": [ "Function([-9.55804991e-01, -5.29095427e-15, -3.05007135e+00,\n", " -4.73748833e-15, 9.51428681e-01, -5.30012041e-15,\n", " -9.13950067e-02, -4.67907612e-15, 4.37386282e-03,\n", " -4.20804838e-15, -1.26261106e-04, -4.37960451e-15,\n", " 2.44435655e-06, -2.91569988e-15, -3.46863840e-08,\n", " 0.00000000e+00])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u = TrialFunction(SN)\n", "v = TestFunction(SN)\n", "D0 = inner(u, v)\n", "D1 = inner(Dx(u, 0, 1), v)\n", "D2 = inner(Dx(u, 0, 2), v)\n", "w0 = Function(SN) # work array\n", "uN1 = Function(SN)\n", "uN2 = Function(SN)\n", "uN1 = D0.solve(D1.matvec(uN, w0), uN1)\n", "uN2 = D0.solve(D2.matvec(uN, w0), uN2)\n", "uN1" ] }, { "cell_type": "markdown", "id": "8203a97b", "metadata": { "editable": true }, "source": [ "We can look at the sparsity patterns of $(d^{(1)}_{kj})$ and $(d^{(2)}_{kj})$" ] }, { "cell_type": "code", "execution_count": 5, "id": "0288b4fe", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.056503Z", "iopub.status.busy": "2024-02-19T13:48:57.056318Z", "iopub.status.idle": "2024-02-19T13:48:57.701357Z", "shell.execute_reply": "2024-02-19T13:48:57.700856Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAESCAYAAABdK7eSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgGUlEQVR4nO3df2xV9f3H8delhcuPtNe1BuodLdSkAy2KCGoUVMi0CyMoLuoErUyTRWMRKoZBp07itFd0Y7h1QMof6kJQ/phUZma027BonAotVcc2kNm1VSQdi7m34LhAe75/EO6XQml723M+58d9PpIbcm/v7etz6fm8827PvfcdsizLEgAAgCHD3F4AAADILDQfAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRNB8AAMAoXzUf69evV3FxsUaOHKnp06fr3XffdSQnFovpqquuUk5OjsaOHasFCxZo3759jmSdLz8UCqmystLRnC+//FL33HOP8vPzNXr0aF1xxRVqbGx0JOvkyZN6/PHHVVxcrFGjRuniiy/WU089pe7ublu+/86dOzV//nxFo1GFQiHV1dX1+LplWVq9erWi0ahGjRql2bNna+/evY7knThxQitXrtRll12mMWPGKBqN6t5779XBgwcdyTvbAw88oFAopHXr1g06L2ioHfYJUt2QqB1nMlk7fNN8bN26VZWVlXrssce0Z88eXX/99Zo7d67a2tpsz2poaFBFRYU++OAD1dfX6+TJkyorK9PRo0dtzzrbrl27VFtbq8svv9zRnK+//lozZ87U8OHD9eabb+rvf/+7fvnLX+qCCy5wJG/NmjXauHGjampq9I9//EPPPfecnn/+ef3mN7+x5fsfPXpUU6dOVU1NTa9ff+6557R27VrV1NRo165dKigo0M0336zOzk7b87755hs1NTXpiSeeUFNTk1577TXt379ft9xyy6Cy+ss7U11dnT788ENFo9FBZwUNtcM+QasbErXjNOO1w/KJq6++2nrwwQd73DZ58mRr1apVjmd3dHRYkqyGhgZHczo7O62SkhKrvr7euvHGG61ly5Y5lrVy5Upr1qxZjn3/s82bN8+6//77e9z2gx/8wLrnnntsz5Jkbdu2LXW9u7vbKigosJ599tnUbceOHbMikYi1ceNG2/N689FHH1mSrNbWVsfyvvjiC+vb3/629be//c2aMGGC9atf/WrIWUFA7bBPkOuGZVE7TNYOX/zl4/jx42psbFRZWVmP28vKyvT+++87nh+PxyVJeXl5juZUVFRo3rx5uummmxzNkaTt27drxowZuuOOOzR27FhNmzZNmzZtcixv1qxZ+vOf/6z9+/dLkj7++GO99957+v73v+9Y5mktLS06dOhQj+MnHA7rxhtvNHL8SKeOoVAo5NhviN3d3SovL9eKFStUWlrqSIYfUTvslUl1Q6J2OCnbWNIQHD58WF1dXRo3blyP28eNG6dDhw45mm1ZlpYvX65Zs2ZpypQpjuW8+uqrampq0q5duxzLONPnn3+uDRs2aPny5frpT3+qjz76SEuXLlU4HNa9995re97KlSsVj8c1efJkZWVlqaurS88884wWLlxoe9bZTh8jvR0/ra2tjucfO3ZMq1at0qJFi5Sbm+tIxpo1a5Sdna2lS5c68v39itphr0yqGxK1w0m+aD5OC4VCPa5blnXObXZbsmSJPvnkE7333nuOZbS3t2vZsmV6++23NXLkSMdyztTd3a0ZM2aourpakjRt2jTt3btXGzZscKSIbN26VZs3b9aWLVtUWlqq5uZmVVZWKhqNavHixbbn9caN4+fEiRO666671N3drfXr1zuS0djYqBdeeEFNTU2OPx+/onbYIxPrhkTtcIIvTrtceOGFysrKOuc3lY6OjnM6Ujs9/PDD2r59u3bs2KHx48c7ltPY2KiOjg5Nnz5d2dnZys7OVkNDg379618rOztbXV1dtmdedNFFuvTSS3vcdskllzjyIjxJWrFihVatWqW77rpLl112mcrLy/XII48oFos5knemgoICSTJ+/Jw4cUJ33nmnWlpaVF9f79hvLu+++646OjpUVFSUOn5aW1v16KOPauLEiY5k+gW1w97akUl1Q6J2OMkXzceIESM0ffp01dfX97i9vr5e1113ne15lmVpyZIleu211/SXv/xFxcXFtmec6bvf/a4+/fRTNTc3py4zZszQ3XffrebmZmVlZdmeOXPmzHPeArh//35NmDDB9izp1Ku4hw3rebhlZWXZ+pa58ykuLlZBQUGP4+f48eNqaGhw5PiR/r94fPbZZ/rTn/6k/Px8R3Ikqby8XJ988kmP4ycajWrFihV66623HMv1A2qHvbUjk+qGRO1wkm9Ouyxfvlzl5eWaMWOGrr32WtXW1qqtrU0PPvig7VkVFRXasmWLXn/9deXk5KS63kgkolGjRtmel5OTc8454TFjxig/P9+xc8WPPPKIrrvuOlVXV+vOO+/URx99pNraWtXW1jqSN3/+fD3zzDMqKipSaWmp9uzZo7Vr1+r++++35fsfOXJEBw4cSF1vaWlRc3Oz8vLyVFRUpMrKSlVXV6ukpEQlJSWqrq7W6NGjtWjRItvzotGobr/9djU1NemNN95QV1dX6hjKy8vTiBEjbH9+Zxeo4cOHq6CgQJMmTRrU8wsSaod9glY3JGqHa7XD8ffT2Oi3v/2tNWHCBGvEiBHWlVde6djb1yT1ennxxRcdyeuN02+1tSzL+sMf/mBNmTLFCofD1uTJk63a2lrHshKJhLVs2TKrqKjIGjlypHXxxRdbjz32mJVMJm35/jt27Oj1Z7Z48WLLsk69Ze7JJ5+0CgoKrHA4bN1www3Wp59+6kheS0vLeY+hHTt2OPL8zsZbbXuidtgnSHXDsqgdZzNVO0KWZVl2NjMAAAB98cVrPgAAQHDQfAAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMIrmAwAAGOW75iOZTGr16tVKJpPkkee5PDcy3XiOfhP0nwnHOXlezzub7z5kLJFIKBKJKB6POzZshzzy/JTpxnP0m6D/TDjOyfN63tl895cPAADgbzQfAADAKM9Nte3u7tbBgweVk5OjUCh0ztcTiUSPf51GHnlez3Qiz7IsdXZ2KhqNnjPS3Kv6qh1B+Jl4Kc+NTPK8n5dO3fDcaz6++OILFRYWur0MAJLa29s1fvx4t5cxINQOwBsGUjc895ePnJwcSacWz4vnAHckEgkVFham9qMfUDsAd6VTNxxrPtavX6/nn39eX331lUpLS7Vu3Tpdf/31/T7u9J9Lc3NzKSCAy3o79emkwdYNidoBeMVA6oYjJ3O3bt2qyspKPfbYY9qzZ4+uv/56zZ07V21tbU7EAQgA6gaQORx5zcc111yjK6+8Uhs2bEjddskll2jBggWKxWJ9Ptbt9x4DcGcfDqVuSNQOwG3p7EHb//Jx/PhxNTY2qqysrMftZWVlev/99+2OG7i2Nqmp6dS/5JHnxcyg5/XBs3VjkDLhRxn050ies2x/zcfhw4fV1dWlcePG9bh93LhxOnTo0Dn3TyaTPT7e1ZG3GbW1SZMmSceOSSNHSvv2SUVF9ueQF4w8NzKDntePdOuGZKh2DEIm/CiD/hzJc55jb+A/+wUnlmX1+iKUWCymSCSSujjyVrnDh0/9L0un/j182P4M8oKT50Zm0PMGaKB1QzJUOwYhE36UQX+O5DnP9ubjwgsvVFZW1jm/rXR0dJzzW40kVVVVKR6Ppy7t7e12L0m68MJT7Z106t8LL7Q/g7zg5LmRGfS8fpeTXt2QDNWOQciEH2XQnyN5zrP9tMuIESM0ffp01dfX67bbbkvdXl9fr1tvvfWc+4fDYYXDYbuX0VNR0am/Kx0+fOp/2em/L5Hn7zw3MoOe149064ZkqHYMQib8KIP+HMlzniPvdtm6davKy8u1ceNGXXvttaqtrdWmTZu0d+9eTZgwoc/H8op1wH1u7MOh1A2J2gG4LZ096MiHjP3whz/Uf//7Xz311FP66quvNGXKFP3xj38cUAEBkJmoG0Dm8NxsF357Adznx33oxzUDQeLq53wAAAD0heYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRjnzOR6C0tZn9GDjy/J3nRqYbzxGO4HAlL1PyaD76EvRRg+T5P9ML4ylhCw5X8jIpj9MufQn6qEHy/J/phfGUsAWHK3mZlEfz0Zegjxokz/+ZXhhPCVtwuJKXSXl8vHp//HpCjTx38tzIdCDPc/twAPy45rNxuJLn57x09iDNB4Bz+HEf+nHNQJAw2wUAAHgWzQcAADCK5gMAABhF8wEAAIyi+QAAAEbZ3nzEYjFdddVVysnJ0dixY7VgwQLt27fP7hgAAULdADKL7c1HQ0ODKioq9MEHH6i+vl4nT55UWVmZjh49ancUgICgbgCZxfHP+fjPf/6jsWPHqqGhQTfccEO/9+e9+oD73N6H6dYNyf01A5kunT3o+GC5eDwuScrLy3M6ylu89LFz5PkjM+h5acjYujFIbA/yvJ53NkebD8uytHz5cs2aNUtTpkzp9T7JZFLJZDJ1PZFIOLkkM/w8apA8dzKDnpeGgdQNKaC1YxDYHuR5Pa83jr7bZcmSJfrkk0/0yiuvnPc+sVhMkUgkdSksLHRySWb4edQgee5kBj0vDQOpG1JAa8cgsD3I83pebxxrPh5++GFt375dO3bs0Pjx4897v6qqKsXj8dSlvb3dqSWZ4+dRg+S5kxn0vAEaaN2QAlo7BoHtQZ7X83pj+wtOLcvSww8/rG3btumdd95RSUlJWo8PzIvGgn4CL+h5bmR6KM/0Phxq3ZACVDsGge1BnhfyXJ1q+9BDD2nLli16/fXXNWnSpNTtkUhEo0aN6vfxmVxAAK8wvQ+HWjckagfgNlebj1Ao1OvtL774on70ox/1+3gKCOA+0/twqHVDonYAbnP1rbYOf2wIgACibgCZhdkuAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMovkAAABGOT7VFmkKwsfcZXKeG5luj6eErwX9cCXPm3k0H14S9NGGQc9zI9ML4ynhW0E/XMnzbh6nXbwk6KMNg57nRqYXxlPCt4J+uJLn3TyaDy8J+mjDoOe5kemF8ZTwraAfruR5N8/22S5DlfHzGfx6Ao889zIdyPPjPvTjmr0gAIcreR7Jc3Ww3FBRQAD3+XEf+nHNQJCkswc57QIAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYJTjzUcsFlMoFFJlZaXTUQACgroBBJujzceuXbtUW1uryy+/3MkYAAFC3QCCz7Hm48iRI7r77ru1adMmfetb33IqBqe1tUlNTaf+Jc+fmUHPGwDqhn8E/XAlz1mOTbWtqKjQvHnzdNNNN+npp58+7/2SyaSSyWTqeiKRcGpJweXn0YZezHMjM+h5AzTQuiFRO9wU9MOVPOc58pePV199VU1NTYrFYv3eNxaLKRKJpC6FhYVOLCnY/Dza0It5bmQGPW8A0qkbErXDTUE/XMlznu3NR3t7u5YtW6bNmzdr5Onxd32oqqpSPB5PXdrb2+1eUvD5ebShF/PcyAx6Xj/SrRsStcNNQT9cyXOe7YPl6urqdNtttykrKyt1W1dXl0KhkIYNG6ZkMtnja2djONQgeWm0YRDy3Mj0UJ7pfTjUuiFRO0zz0OFKnkfyXJ1q29nZqdbW1h633XfffZo8ebJWrlypKVOm9Pl4CgjgPtP7cKh1Q6J2AG5LZw/a/oLTnJyccwrFmDFjlJ+fP6ACAiDzUDeAzMInnAIAAKMce6vtmd555x0TMQAChLoBBBd/+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYJSRd7vAw4LwsXpeynMj043nCAxS0LcHeQND85HJgj5KkWm4gKcEfXuQN3CcdslkQR+lyDRcwFOCvj3IGziaj0wW9FGKTMMFPCXo24O8gbN9sNxQMRzKML+eMPRqnhuZDuT5cR/6cc2ZKADbg7zz5Lk61XaoKCCA+/y4D/24ZiBI0tmDnHYBAABG0XwAAACjaD4AAIBRNB8AAMAomg8AAGAUzQcAADDKkebjyy+/1D333KP8/HyNHj1aV1xxhRobG52IAhAQ1A0gc9g+2+Xrr7/WzJkzNWfOHL355psaO3as/vWvf+mCCy6wOwpAQFA3gMxie/OxZs0aFRYW6sUXX0zdNnHiRLtj4DYvfaxeUDKDntcH6gb6E/TtEfS8s9l+2mX79u2aMWOG7rjjDo0dO1bTpk3Tpk2bznv/ZDKpRCLR4wKPOz3acPr0U/+2tQUrz43MoOf1I926IVE7MknQt0fQ83pje/Px+eefa8OGDSopKdFbb72lBx98UEuXLtXvfve7Xu8fi8UUiURSl8LCQruXBLv5eZSiVzODntePdOuGRO3IJEHfHkHP643tzUd3d7euvPJKVVdXa9q0aXrggQf04x//WBs2bOj1/lVVVYrH46lLe3u73UuC3fw8StGrmUHP60e6dUOidmSSoG+PoOf1xvbXfFx00UW69NJLe9x2ySWX6Pe//32v9w+HwwqHw3YvA04qKpL27TN3wtB0nhuZQc/rR7p1Q6J2ZJKgb4+g5/XG9uZj5syZ2rdvX4/b9u/frwkTJtgdBTcVFZk9Yk3nuZEZ9Lw+UDfQn6Bvj6Dnnc320y6PPPKIPvjgA1VXV+vAgQPasmWLamtrVVFRYXcUgICgbgCZxfbm46qrrtK2bdv0yiuvaMqUKfr5z3+udevW6e6777Y7CkBAUDeAzBKyLMtyexFnSiQSikQiisfjys3NdXs5QEby4z7045qBIElnDzLbBQAAGEXzAQAAjKL5AAAARtF8AAAAo2g+AACAUbZ/yBjQp0wY3ZgJzxHwiaBvR7/m0XzAnNOjFI8dOzVQYN8+Z3eL6Tw3Mt14joBPBH07+jmP0y4wJxNGN2bCcwR8Iujb0c95NB8wJxNGN2bCcwR8Iujb0c95nHaBOZkwujETniPgE0Hfjn7O4+PVAZzDj/vQj2sGgoSPVwcAAJ5F8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCjbm4+TJ0/q8ccfV3FxsUaNGqWLL75YTz31lLq7u+2OAhAQ1A0gs9j+IWNr1qzRxo0b9fLLL6u0tFS7d+/Wfffdp0gkomXLltkdByAAqBtAZrG9+fjrX/+qW2+9VfPmzZMkTZw4Ua+88op2795tdxQyCdNp/Z/XB+oGvCbo29Ht7W978zFr1ixt3LhR+/fv13e+8x19/PHHeu+997Ru3bpe759MJpVMJlPXE4mE3UuC3zGd1v95/Ui3bkjUDjgn6NvRC9vf9td8rFy5UgsXLtTkyZM1fPhwTZs2TZWVlVq4cGGv94/FYopEIqlLYWGh3UuC3zGd1v95/Ui3bkjUDjgn6NvRC9vf9uZj69at2rx5s7Zs2aKmpia9/PLL+sUvfqGXX3651/tXVVUpHo+nLu3t7XYvCX7HdFr/5/Uj3bohUTvgnKBvRy9sf9sHyxUWFmrVqlWqqKhI3fb0009r8+bN+uc//9nv4xkOhV7xmg+jeab34VDrhkTtgL08tB19k5fOHrT9NR/ffPONhg3r+QeVrKws3jKHoSkqMn9S0nRm0PP6QN2A1wR9O7q9/W1vPubPn69nnnlGRUVFKi0t1Z49e7R27Vrdf//9dkcBCAjqBpBZbD/t0tnZqSeeeELbtm1TR0eHotGoFi5cqJ/97GcaMWJEv4/nT6eA+0zvw6HWDYnaAbgtnT1oe/MxVBQQwH1+3Id+XDMQJOnsQWa7AAAAo2g+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAYZfuHjAGewseyAzAo6NvfrjyaDwSXG3OjM3E2NgBJwd/+duZx2gXB5cbc6EycjQ1AUvC3v515NB8ILjfmRmfibGwAkoK//e3M47QLgquo6NTfBU2eEDWd6cZzBNCroG9/O/OY7QLgHH7ch35cMxAkzHYBAACeRfMBAACMovkAAABG0XwAAACjaD4AAIBRaTcfO3fu1Pz58xWNRhUKhVRXV9fj65ZlafXq1YpGoxo1apRmz56tvXv32rVeAD5E3QBwprSbj6NHj2rq1Kmqqanp9evPPfec1q5dq5qaGu3atUsFBQW6+eab1dnZOeTFAvAn6gaAHqwhkGRt27Ytdb27u9sqKCiwnn322dRtx44dsyKRiLVx48YBfc94PG5JsuLx+FCWBmAInNyHTtQNy6J2AG5LZw/a+pqPlpYWHTp0SGVlZanbwuGwbrzxRr3//vt2RgHOa2uTmppO/UueY6gbQOZtf1s/Xv3QoUOSpHHjxvW4fdy4cWptbe31MclkUslkMnU9kUjYuSRgcPw8LtKLeX0YTN2QqB0Ijkzc/o682yUUCvW4blnWObedFovFFIlEUpfCwkInlgSkx8/jIr2YNwDp1A2J2oHgyMTtb2vzUVBQIOn/f5M5raOj45zfak6rqqpSPB5PXdrb2+1cEjA4fh4X6cW8PgymbkjUDgRHJm5/W0+7FBcXq6CgQPX19Zo2bZok6fjx42poaNCaNWt6fUw4HFY4HLZzGcDQ+XlcpBfz+jCYuiFROxAcmbj9024+jhw5ogMHDqSut7S0qLm5WXl5eSoqKlJlZaWqq6tVUlKikpISVVdXa/To0Vq0aJGtCwccV1RkdlcGOI+6AfQtwNu/V2k3H7t379acOXNS15cvXy5JWrx4sV566SX95Cc/0f/+9z899NBD+vrrr3XNNdfo7bffVk5Ojn2rBuAr1A0AZwpZlmW5vYgzJRIJRSIRxeNx5ebmur0cICP5cR/6cc1AkKSzB5ntAgAAjKL5AAAARtF8AAAAo2g+AACAUTQfAADAKJoPAABglK2fcApAp6Y2mfzoQNN5ADzDr+WG5gOwUyaOpwTgCj+XG067AHbKxPGUAFzh53JD8wHYKRPHUwJwhZ/LDaddADtl4nhKAK7wc7mh+QDslmnjKQG4xq/lhtMuAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMSrv52Llzp+bPn69oNKpQKKS6urrU106cOKGVK1fqsssu05gxYxSNRnXvvffq4MGDdq4ZgM9QNwCcKe3m4+jRo5o6dapqamrO+do333yjpqYmPfHEE2pqatJrr72m/fv365ZbbrFlsQD8iboB4Ewhy7KsQT84FNK2bdu0YMGC895n165duvrqq9Xa2qqiAbw5OJFIKBKJKB6PKzc3d7BLAzAETu5DJ+qGRO0A3JbOHnT8Q8bi8bhCoZAuuOACp6MAf/PreEoHUDcAZ7m9/R1tPo4dO6ZVq1Zp0aJF5+2Cksmkkslk6noikXBySYA3+Xk8pc0GUjckagcwWF7Y/o692+XEiRO666671N3drfXr15/3frFYTJFIJHUpLCx0akmAd/l5PKWNBlo3JGoHMFhe2P6ONB8nTpzQnXfeqZaWFtXX1/f520tVVZXi8Xjq0t7e7sSSAG/z83hKm6RTNyRqBzBYXtj+tp92OV1APvvsM+3YsUP5+fl93j8cDiscDtu9DMBf/Dye0gbp1g2J2gEMlhe2f9rNx5EjR3TgwIHU9ZaWFjU3NysvL0/RaFS33367mpqa9MYbb6irq0uHDh2SJOXl5WnEiBH2rRwIGr+OpxwA6gbgLW4Pw077rbbvvPOO5syZc87tixcv1urVq1VcXNzr43bs2KHZs2f3+/15uxzgPrv3odN1Q6J2AG5z9K22s2fPVl/9yhA+NgRAQFE3AJyJ2S4AAMAomg8AAGAUzQcAADCK5gMAABhF8wEAAIyi+QAAAEY5PtUWgMPcHk8JIGPYVW5oPgA/88J4SgAZwc5yw2kXwM+8MJ4SQEaws9zQfAB+5oXxlAAygp3lhtMugJ95YTwlgIxgZ7mh+QD8zu3xlAAyhl3lhtMuAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMovkAAABGpd187Ny5U/Pnz1c0GlUoFFJdXd157/vAAw8oFApp3bp1Q1giAL+jbgA4U9rNx9GjRzV16lTV1NT0eb+6ujp9+OGHikajg14cgGCgbgA4U9ofMjZ37lzNnTu3z/t8+eWXWrJkid566y3Nmzdv0IsD4CCD03CpGwDOZPsnnHZ3d6u8vFwrVqxQaWlpv/dPJpNKJpOp64lEwu4lATibx6bhpls3JGoH4Ge2v+B0zZo1ys7O1tKlSwd0/1gspkgkkroUFhbavSQAZ/PYNNx064ZE7QD8zNbmo7GxUS+88IJeeuklhUKhAT2mqqpK8Xg8dWlvb7dzSQB646FpuIOpGxK1A/AzW5uPd999Vx0dHSoqKlJ2drays7PV2tqqRx99VBMnTuz1MeFwWLm5uT0uABx2ejxlY6Prp1wGUzckagfgZ7a+5qO8vFw33XRTj9u+973vqby8XPfdd5+dUQCGyiPTcKkbQOZJu/k4cuSIDhw4kLre0tKi5uZm5eXlqaioSPn5+T3uP3z4cBUUFGjSpElDXy0AX6JuADhT2s3H7t27NWfOnNT15cuXS5IWL16sl156ybaFAQgO6gaAM6XdfMyePVuWZQ34/v/+97/TjQAQMNQNAGditgsAADCK5gMAABhF8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMIrmAwAAGEXzAQAAjKL5AAAARtF8AAAAo9Keauu005MvE4mEyysBMtfp/ZfOJFq3UTsAd6VTNzzXfHR2dkqSCgsLXV4JgM7OTkUiEbeXMSDUDsAbBlI3QpbHfrXp7u7WwYMHlZOTo1AodM7XE4mECgsL1d7ertzcXMfXQx55Xs90Is+yLHV2dioajWrYMH+cne2rdgThZ+KlPDcyyfN+Xjp1w3N/+Rg2bJjGjx/f7/1yc3ONbTLyyPNDpt15fvmLx2kDqR1+/5l4Lc+NTPK8nTfQuuGPX2kAAEBg0HwAAACjfNd8hMNhPfnkkwqHw+SR57k8NzLdeI5+E/SfCcc5eV7PO5vnXnAKAACCzXd/+QAAAP5G8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMOr/AFgVLw8t6dWaAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "fig, (ax1, ax2) = plt.subplots(1, 2)\n", "ax1.spy(D1.diags(), markersize=2, color='r')\n", "ax2.spy(D2.diags(), markersize=2, color='b')" ] }, { "cell_type": "markdown", "id": "798dcfdd", "metadata": { "editable": true }, "source": [ "just to see that they are upper triangular. We now ask is there a better and faster\n", "way to get `uN1` and `uN2`? A better approach would involve only sparse\n", "matrices, like the diagonal $(d^{(0)}_{kj})$. But how do we get there?\n", "Most textbooks on spectral methods use fast recursive methods to\n", "find the coefficients $\\{\\hat{u}_k^{(n)}\\}$. Here we will show a fast Galerkin approach.\n", "\n", "It turns out that a simple change of test space/function will be sufficient.\n", "Let us first replace the test space $S_N$ with\n", "the Dirichlet space $D_N=\\{v \\in S_N | v(\\pm 1) = 0\\}$ using basis\n", "functions $v=T_k-T_{k+2}$ and see\n", "what happens. Because of the two boundary conditions,\n", "the number of degrees of freedom is reduced by two, and we need to use a\n", "space with $N+2$ quadrature points in order to get a square matrix system.\n", "The method now becomes classified as Chebyshev-Petrov-Galerkin, as we\n", "wish to find $u_N^{(1)} \\in S_N$ such that" ] }, { "cell_type": "markdown", "id": "a0edfd85", "metadata": { "editable": true }, "source": [ "$$\n", "(u_N^{(1)}-u^{(1)}, v)_{\\omega^{-1/2}} = 0, \\quad \\forall v \\in D_{N+2}.\n", "$$" ] }, { "cell_type": "markdown", "id": "6fb9034a", "metadata": { "editable": true }, "source": [ "The implementation is straightforward" ] }, { "cell_type": "code", "execution_count": 6, "id": "78eb36cf", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.704227Z", "iopub.status.busy": "2024-02-19T13:48:57.703893Z", "iopub.status.idle": "2024-02-19T13:48:57.721596Z", "shell.execute_reply": "2024-02-19T13:48:57.720958Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.00000000e+00 -7.88860905e-31 0.00000000e+00 0.00000000e+00\n", " 1.11022302e-16 -7.88860905e-31 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 4.23516474e-22 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n" ] } ], "source": [ "DN = FunctionSpace(N+2, 'C', bc=(0, 0))\n", "v = TestFunction(DN)\n", "D0 = inner(u, v)\n", "D1 = inner(Dx(u, 0, 1), v)\n", "uN11 = Function(SN)\n", "uN11 = D0.solve(D1.matvec(uN, w0), uN11)\n", "print(uN11-uN1)" ] }, { "cell_type": "markdown", "id": "6fc6f39f", "metadata": { "editable": true }, "source": [ "and since `uN11 = uN1` we see that we have achived the same result as in\n", "the regular projection. However, the matrices in use now look like" ] }, { "cell_type": "code", "execution_count": 7, "id": "d5e3bb84", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:57.724288Z", "iopub.status.busy": "2024-02-19T13:48:57.724086Z", "iopub.status.idle": "2024-02-19T13:48:58.016969Z", "shell.execute_reply": "2024-02-19T13:48:58.016477Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAESCAYAAABdK7eSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAexUlEQVR4nO3df2zU9R3H8ddR7PEj7bmWQL1xxZoQ0aKoVI2AClG7IMGxRZ2ilek/khWh1jDo1Mnc6IluDLcOSP1DWAjKH5PKzIx2DosGf5SWqnMbyOxoFZuOxdwVXA9ov/uj6Y3S0vba7/dz3+/d85Fcmvv22ven6X1eed997+7tsyzLEgAAgCFjkr0AAACQXmg+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRnmo+Nm/erIKCAo0bN06zZ8/WO++840idcDisa6+9VllZWZo8ebKWLFmiQ4cOOVLrfPV9Pp/KysocrfPll1/q/vvvV25uriZMmKCrrrpKDQ0NjtQ6c+aMnnjiCRUUFGj8+PG65JJL9PTTT6u7u9uW379v3z4tXrxYwWBQPp9PNTU1fb5vWZbWrVunYDCo8ePHa/78+fr0008dqXf69GmtWbNGV1xxhSZOnKhgMKgHHnhAx44dc6TeuR5++GH5fD5t2rRpxPVSDdlhn1TKDYnsOJvJ7PBM87Fr1y6VlZXp8ccf18GDB3XjjTdq4cKFamlpsb1WXV2dSktL9f7776u2tlZnzpxRcXGxTp48aXutc9XX16u6ulpXXnmlo3W+/vprzZ07VxdccIFef/11/e1vf9OvfvUrXXjhhY7U27Bhg7Zu3aqqqir9/e9/17PPPqvnnntOv/3tb235/SdPntSsWbNUVVU14PefffZZbdy4UVVVVaqvr1deXp5uu+02dXR02F7vm2++UWNjo5588kk1NjbqlVde0eHDh3XHHXeMqNZQ9c5WU1OjDz74QMFgcMS1Ug3ZYZ9Uyw2J7OhlPDssj7juuuus5cuX9zk2Y8YMa+3atY7Xbm9vtyRZdXV1jtbp6Oiwpk+fbtXW1lo333yztWrVKsdqrVmzxpo3b55jv/9cixYtsh566KE+x77//e9b999/v+21JFm7d++OX+/u7rby8vKsZ555Jn6ss7PTCgQC1tatW22vN5APP/zQkmQdPXrUsXpffPGF9e1vf9v661//ak2bNs369a9/PepaqYDssE8q54ZlkR0ms8MTz3ycOnVKDQ0NKi4u7nO8uLhY+/fvd7x+JBKRJOXk5Dhap7S0VIsWLdKtt97qaB1J2rNnj4qKinTXXXdp8uTJuvrqq/XCCy84Vm/evHl66623dPjwYUnSRx99pHfffVe33367YzV7NTc3q62trc/9x+/36+abbzZy/5F67kM+n8+xR4jd3d0qKSnR6tWrVVhY6EgNLyI77JVOuSGRHU4aa6zSKBw/flxdXV2aMmVKn+NTpkxRW1ubo7Uty1J5ebnmzZunmTNnOlbn5ZdfVmNjo+rr6x2rcbbPP/9cW7ZsUXl5uX7yk5/oww8/1MqVK+X3+/XAAw/YXm/NmjWKRCKaMWOGMjIy1NXVpfXr1+vee++1vda5eu8jA91/jh496nj9zs5OrV27VkuXLlV2drYjNTZs2KCxY8dq5cqVjvx+ryI77JVOuSGRHU7yRPPRy+fz9bluWVa/Y3ZbsWKFPv74Y7377ruO1WhtbdWqVav05ptvaty4cY7VOVt3d7eKiopUWVkpSbr66qv16aefasuWLY6EyK5du7Rjxw7t3LlThYWFampqUllZmYLBoJYtW2Z7vYEk4/5z+vRp3XPPPeru7tbmzZsdqdHQ0KDnn39ejY2Njv89XkV22CMdc0MiO5zgidMukyZNUkZGRr9HKu3t7f06Ujs98sgj2rNnj/bu3aupU6c6VqehoUHt7e2aPXu2xo4dq7Fjx6qurk6/+c1vNHbsWHV1ddle86KLLtLll1/e59hll13myIvwJGn16tVau3at7rnnHl1xxRUqKSnRo48+qnA47Ei9s+Xl5UmS8fvP6dOndffdd6u5uVm1tbWOPXJ555131N7ervz8/Pj95+jRo3rsscd08cUXO1LTK8gOe7MjnXJDIjuc5InmIzMzU7Nnz1ZtbW2f47W1tZozZ47t9SzL0ooVK/TKK6/oL3/5iwoKCmyvcbZbbrlFn3zyiZqamuKXoqIi3XfffWpqalJGRobtNefOndvvLYCHDx/WtGnTbK8l9byKe8yYvne3jIwMW98ydz4FBQXKy8vrc/85deqU6urqHLn/SP8Pj88++0x//vOflZub60gdSSopKdHHH3/c5/4TDAa1evVqvfHGG47V9QKyw97sSKfckMgOJ3nmtEt5eblKSkpUVFSkG264QdXV1WppadHy5cttr1VaWqqdO3fq1VdfVVZWVrzrDQQCGj9+vO31srKy+p0TnjhxonJzcx07V/zoo49qzpw5qqys1N13360PP/xQ1dXVqq6udqTe4sWLtX79euXn56uwsFAHDx7Uxo0b9dBDD9ny+0+cOKEjR47Erzc3N6upqUk5OTnKz89XWVmZKisrNX36dE2fPl2VlZWaMGGCli5danu9YDCoO++8U42NjXrttdfU1dUVvw/l5OQoMzPT9r/v3IC64IILlJeXp0svvXREf18qITvsk2q5IZEdScsOx99PY6Pf/e531rRp06zMzEzrmmuucezta5IGvLz44ouO1BuI02+1tSzL+uMf/2jNnDnT8vv91owZM6zq6mrHakWjUWvVqlVWfn6+NW7cOOuSSy6xHn/8cSsWi9ny+/fu3Tvg/2zZsmWWZfW8Ze6pp56y8vLyLL/fb910003WJ5984ki95ubm896H9u7d68jfdy7eatsX2WGfVMoNyyI7zmUqO3yWZVl2NjMAAACD8cRrPgAAQOqg+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMMpzzUcsFtO6desUi8WoRz3X1UtGzWT8jV6T6v8T7ufUc3u9c3nuQ8ai0agCgYAikYhjw3aoRz0v1UzG3+g1qf4/4X5OPbfXO5fnnvkAAADeRvMBAACMct1U2+7ubh07dkxZWVny+Xz9vh+NRvt8dRr1qOf2mk7UsyxLHR0dCgaD/Uaau9Vg2ZEK/xM31UtGTeq5v14iueG613x88cUXCoVCyV4GAEmtra2aOnVqspcxLGQH4A7DyQ3XPfORlZUlqWfxvHgOSI5oNKpQKBTfj15AdgDJlUhuONZ8bN68Wc8995y++uorFRYWatOmTbrxxhuH/Lnep0uzs7MJECDJBjr16aSR5oZEdgBuMZzccORk7q5du1RWVqbHH39cBw8e1I033qiFCxeqpaXFiXIAUgC5AaQPR17zcf311+uaa67Rli1b4scuu+wyLVmyROFweNCfTfZ7jwEkZx+OJjcksgNItkT2oO3PfJw6dUoNDQ0qLi7uc7y4uFj79+8ffYGWFqmxseerCabrAWnI8dwYIbY/4AzbX/Nx/PhxdXV1acqUKX2OT5kyRW1tbf1uH4vF+ny866Bv+2lpkS69VOrslMaNkw4dkvLzbVt70usBaSrR3JASzI4RYPsDznHsDfznvuDEsqwBX4QSDocVCATil0HfKnf8eE8SSD1fjx+3c8nJrwekueHmhpRgdowA2x9wju3Nx6RJk5SRkdHv0Up7e3u/RzWSVFFRoUgkEr+0trYO9st7HoJIPV8nTbJz6cmvB6SpRHNDSjA7RrQmtj/gFNubj8zMTM2ePVu1tbV9jtfW1mrOnDn9bu/3++NvjRvyLXL5+T3PfTY0mHkO1HQ9IE0lmhtSgtkxAmx/wDmOfM5HeXm5SkpKVFRUpBtuuEHV1dVqaWnR8uXLR//L8/PNpoDpekCacjQ3RojtDzjDkebjBz/4gf7zn//o6aef1ldffaWZM2fqT3/6k6ZNm+ZEOQApgNwA0ofrZrvwXn0g+by4D724ZiCVJPVzPgAAAAZD8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwChHPufDlVpaeoYzTJpk5lODTNcD4Bpsf2Bw6dF8MA0XgCFsf2Bo6XHahWm4AAxh+wNDS4/mg2m4AAxh+wNDS4/TLr3jKU2dhDVdD4BrsP2BoaVH8yExDReAMWx/YHDpcdoFAAC4Bs0HAAAwiuYDAAAYRfMBAACMovkAAABG2d58hMNhXXvttcrKytLkyZO1ZMkSHTp0yO4yAFIIuQGkF9ubj7q6OpWWlur9999XbW2tzpw5o+LiYp08edLuUgBSBLkBpBefZVmWkwX+/e9/a/Lkyaqrq9NNN9005O2j0agCgYAikYiys7OdXBqA80j2Pkw0N6TkrxlId4nsQcc/ZCwSiUiScnJynC7lDKbhAsZ5PjdGiO2PdOFo82FZlsrLyzVv3jzNnDlzwNvEYjHFYrH49Wg06uSSEsM0XMC44eSG5PLsGAG2P9KJo+92WbFihT7++GO99NJL571NOBxWIBCIX0KhkJNLSgzTcAHjhpMbksuzYwTY/kgnjjUfjzzyiPbs2aO9e/dq6tSp571dRUWFIpFI/NLa2urUkhLHNFzAqOHmhuTy7BgBtj/Sie2nXSzL0iOPPKLdu3fr7bffVkFBwaC39/v98vv9di/DHkzDBYxINDckl2fHCLD9kU5sbz5KS0u1c+dOvfrqq8rKylJbW5skKRAIaPz48XaXcx7TcAHHpVxujBDbH+nC9rfa+ny+AY+/+OKL+uEPfzjkz/N2OSD5TO/D0eaGRHYAyZbUt9o6/LEhAFIQuQGkF2a7AAAAo2g+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAY5fhU27TFNFwAhrD94TU0H05gGi4AQ9j+8CJOuziBabgADGH7w4toPpzANFwAhrD94UWcdnEC03ABGML2hxfRfDiFabgADGH7w2s47QIAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYJTjzUc4HJbP51NZWZnTpQCkCHIDSG2ONh/19fWqrq7WlVde6WQZACmE3ABSn2PNx4kTJ3TffffphRde0Le+9S2nyqSelhapsbHnayrWAwZBbphHBCAZHGs+SktLtWjRIt16662D3i4Wiykajfa5pK3e8ZSzZ/d8dToNTNcDhjDc3JDIDjsQAUgWR5qPl19+WY2NjQqHw0PeNhwOKxAIxC+hUMiJJXkD03CRxhLJDYnssAMRgGSxvflobW3VqlWrtGPHDo3rHbU4iIqKCkUikfiltbXV7iV5B9NwkaYSzQ2J7LADEYBk8VmWZdn5C2tqavS9731PGRkZ8WNdXV3y+XwaM2aMYrFYn++dKxqNKhAIKBKJKDs7286leUNLi9nxlKbrwRNM78PR5oZEdowUEQC7JLIHbZ9qe8stt+iTTz7pc+zBBx/UjBkztGbNmiEDJO0xDRdpiNxIHiIAyWB785GVlaWZM2f2OTZx4kTl5ub2Ow4AErkBpBs+4RQAABhl+zMfA3n77bdNlAGQQsgNIHXxzAcAADCK5gMAABhF8wEAAIyi+QAAAEbRfAAAAKOMvNsFBvDJqAAMIgIwGjQfqaB3NGVnZ8+AhkOHnE0D0/UAuAoRgNHitEsqYBouAIOIAIwWzUcqYBouAIOIAIwWp11SQX5+z/Oepk7Amq4HwFWIAIwWzUeqYBouAIOIAIwGp10AAIBRNB8AAMAomg8AAGAUzQcAADCK5gMAABhF8wEAAIxypPn48ssvdf/99ys3N1cTJkzQVVddpYaGBidKAUgR5AaQPmz/nI+vv/5ac+fO1YIFC/T6669r8uTJ+uc//6kLL7zQ7lIAUgS5AaQX25uPDRs2KBQK6cUXX4wfu/jii+0uA7swDRcuQG6kDyIAkgOnXfbs2aOioiLdddddmjx5sq6++mq98MIL5719LBZTNBrtc4EhvaMpZ8/u+drSklr14BmJ5oZEdngREYBetjcfn3/+ubZs2aLp06frjTfe0PLly7Vy5Ur9/ve/H/D24XBYgUAgfgmFQnYvCefDNFy4RKK5IZEdXkQEoJfPsizLzl+YmZmpoqIi7d+/P35s5cqVqq+v13vvvdfv9rFYTLFYLH49Go0qFAopEokoOzvbzqXhXL0PQzo7e0ZTHjrk7POgputhxKLRqAKBgLF9mGhuSGSHFxEBqS2R3LD9NR8XXXSRLr/88j7HLrvsMv3hD38Y8PZ+v19+v9/uZWA4mIYLl0g0NySyw4uIAPSyvfmYO3euDh061OfY4cOHNW3aNLtLwQ5Mw4ULkBvpgwiA5MBrPh599FG9//77qqys1JEjR7Rz505VV1ertLTU7lIAUgS5AaQX25uPa6+9Vrt379ZLL72kmTNn6uc//7k2bdqk++67z+5SAFIEuQGkF9tfcDpapl/oBqA/L+5DL64ZSCWJ7EFmuwAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMMr2DxlDmmAaLgCDiIDUQvOBxDETBoBBREDq4bQLEsc0XAAGEQGph+YDiZs0qefhh9TzddKk1KoHwFWIgNTDaRckjmm4AAwiAlIPzQdGhmm4AAwiAlILp10AAIBRNB8AAMAomg8AAGAUzQcAADCK5gMAABhF8wEAAIyyvfk4c+aMnnjiCRUUFGj8+PG65JJL9PTTT6u7u9vuUgBSBLkBpBfbP+djw4YN2rp1q7Zv367CwkIdOHBADz74oAKBgFatWmV3OQApgNwA0ovtzcd7772n7373u1q0aJEk6eKLL9ZLL72kAwcO2F0KXsQ0XAyA3IBTiAB3sv20y7x58/TWW2/p8OHDkqSPPvpI7777rm6//fYBbx+LxRSNRvtckKJ6R1POnt3ztaUltephxBLNDYnswNCIAPey/ZmPNWvWKBKJaMaMGcrIyFBXV5fWr1+ve++9d8Dbh8Nh/exnP7N7GXCjgUZTOvlQxHQ9jFiiuSGRHRgaEeBetj/zsWvXLu3YsUM7d+5UY2Ojtm/frl/+8pfavn37gLevqKhQJBKJX1pbW+1eEtyCabg4j0RzQyI7MDQiwL18lmVZdv7CUCiktWvXqrS0NH7sF7/4hXbs2KF//OMfQ/58NBpVIBBQJBJRdna2nUuDG/CaD08wvQ9HmxsS2YGBEQHmJLIHbT/t8s0332jMmL5PqGRkZPCWOfRgGi4GQG7AKUSAO9nefCxevFjr169Xfn6+CgsLdfDgQW3cuFEPPfSQ3aUApAhyA0gvtp926ejo0JNPPqndu3ervb1dwWBQ9957r376058qMzNzyJ/nqVMg+Uzvw9HmhkR2AMmWyB60vfkYLQIESD4v7kMvrhlIJYnsQWa7AAAAo2g+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAYZfuHjAGOSMZnJPO5zEDaYvs7i+YD7tc7F7uzs2c61KFDzqdBMmoCcAW2v/M47QL3G2gudirWBOAKbH/n0XzA/ZIxF5tZ3EDaYvs7j9MucL/8/J7nPU2egE1GTQCuwPZ3Hs0HvCEZc7GZxQ2kLba/szjtAgAAjKL5AAAARtF8AAAAo2g+AACAUTQfAADAqISbj3379mnx4sUKBoPy+Xyqqanp833LsrRu3ToFg0GNHz9e8+fP16effmrXegF4ELkB4GwJNx8nT57UrFmzVFVVNeD3n332WW3cuFFVVVWqr69XXl6ebrvtNnV0dIx6sQC8idwA0Ic1CpKs3bt3x693d3dbeXl51jPPPBM/1tnZaQUCAWvr1q3D+p2RSMSSZEUikdEsDcAoOLkPncgNyyI7gGRLZA/a+pqP5uZmtbW1qbi4OH7M7/fr5ptv1v79++0sBQxPS4vU2NjzNZVrehi5gVTC9h8eWz/htK2tTZI0ZcqUPsenTJmio0ePDvgzsVhMsVgsfj0ajdq5JKQzpuF6wkhyQyI74D5s/+Fz5N0uPp+vz3XLsvod6xUOhxUIBOKXUCjkxJKQjpiG6ymJ5IZEdsB92P7DZ2vzkZeXJ+n/j2R6tbe393tU06uiokKRSCR+aW1ttXNJSGdMw/WEkeSGRHbAfdj+w2dr81FQUKC8vDzV1tbGj506dUp1dXWaM2fOgD/j9/uVnZ3d5wLYonc0ZUODuec/k1HT40aSGxLZAfdh+w9fwq/5OHHihI4cORK/3tzcrKamJuXk5Cg/P19lZWWqrKzU9OnTNX36dFVWVmrChAlaunSprQsHhoVpuK5AbiBdsP2HJ+Hm48CBA1qwYEH8enl5uSRp2bJl2rZtm3784x/rv//9r370ox/p66+/1vXXX68333xTWVlZ9q0agKeQGwDO5rMsy0r2Is4WjUYVCAQUiUR4GhVIEi/uQy+uGUgliexBZrsAAACjaD4AAIBRNB8AAMAomg8AAGAUzQcAADCK5gMAABhl62A5IKW0tPQMZ5g0ycynBpmuB8A10m3703wAAzE9npJxmEDaSsftz2kXYCCmx1MyDhNIW+m4/Wk+gIGYHk/JOEwgbaXj9ue0CzCQ3vGUpk7Cmq4HwDXScfvTfADnY3o8JeMwgbSVbtuf0y4AAMAomg8AAGAUzQcAADCK5gMAABhF8wEAAIxKuPnYt2+fFi9erGAwKJ/Pp5qamvj3Tp8+rTVr1uiKK67QxIkTFQwG9cADD+jYsWN2rhmAx5AbAM6WcPNx8uRJzZo1S1VVVf2+980336ixsVFPPvmkGhsb9corr+jw4cO64447bFksAG8iNwCczWdZljXiH/b5tHv3bi1ZsuS8t6mvr9d1112no0ePKn8Yb2KORqMKBAKKRCLKzs4e6dIAjIKT+9CJ3JDIDiDZEtmDjn/IWCQSkc/n04UXXuh0KcAdmIY7auQGMDxe3f6ONh+dnZ1au3atli5det4uKBaLKRaLxa9Ho1EnlwQ4i2m4ozac3JDIDsDL29+xd7ucPn1a99xzj7q7u7V58+bz3i4cDisQCMQvoVDIqSUBzmMa7qgMNzcksgPw8vZ3pPk4ffq07r77bjU3N6u2tnbQRy8VFRWKRCLxS2trqxNLAsxgGu6IJZIbEtkBeHn7237apTdAPvvsM+3du1e5ubmD3t7v98vv99u9DCA5mIY7IonmhkR2AF7e/gk3HydOnNCRI0fi15ubm9XU1KScnBwFg0Hdeeedamxs1Guvvaauri61tbVJknJycpSZmWnfygG3YhpuP+QG4AwPbP8BJfxW27ffflsLFizod3zZsmVat26dCgoKBvy5vXv3av78+UP+ft4uBySf3fvQ6dyQyA4g2Rx9q+38+fM1WL8yio8NAZCiyA0AZ2O2CwAAMIrmAwAAGEXzAQAAjKL5AAAARtF8AAAAo2g+AACAUY5PtQUwTEzDBWBIsrc/zQfgBkzDBWCIG7Y/p10AN2AaLgBD3LD9aT4AN2AaLgBD3LD9Oe0CuAHTcAEY4obtT/MBuAXTcAEYkuztz2kXAABgFM0HAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjEm4+9u3bp8WLFysYDMrn86mmpua8t3344Yfl8/m0adOmUSwRgNeRGwDOlnDzcfLkSc2aNUtVVVWD3q6mpkYffPCBgsHgiBcHIDWQGwDOlvCHjC1cuFALFy4c9DZffvmlVqxYoTfeeEOLFi0a8eIADEOyx1MOA7kBpAa74sb2Tzjt7u5WSUmJVq9ercLCwiFvH4vFFIvF4tej0ajdSwJSlxvGU9og0dyQyA7ANDvjxvYXnG7YsEFjx47VypUrh3X7cDisQCAQv4RCIbuXBKQuN4yntEGiuSGRHYBpdsaNrc1HQ0ODnn/+eW3btk0+n29YP1NRUaFIJBK/tLa22rkkILW5YTzlKI0kNySyAzDNzrix9bTLO++8o/b2duWf9TxMV1eXHnvsMW3atEn/+te/+v2M3++X3++3cxlA+nDDeMpRGkluSGQHYJqdcWNr81FSUqJbb721z7HvfOc7Kikp0YMPPmhnKQC9kj2ecpTIDcA77IqbhJuPEydO6MiRI/Hrzc3NampqUk5OjvLz85Wbm9vn9hdccIHy8vJ06aWXjn61ADyJ3ABwtoSbjwMHDmjBggXx6+Xl5ZKkZcuWadu2bbYtDEDqIDcAnC3h5mP+/PmyLGvYtz/f+VoA6YPcAHA2ZrsAAACjaD4AAIBRNB8AAMAomg8AAGAUzQcAADDK9sFyADzCA9NwAaQmmg8gHaXINFwA3sRpFyAdpcg0XADeRPMBpKMUmIYLwLs47QKkoxSYhgvAu2g+gHTl8Wm4ALyL0y4AAMAomg8AAGCU60679E6+jEajSV4JkL56918ik2iTjewAkiuR3HBd89HR0SFJCoVCSV4JgI6ODgUCgWQvY1jIDsAdhpMbPstlD226u7t17NgxZWVlyefz9ft+NBpVKBRSa2ursrOzHV8P9ajn9ppO1LMsSx0dHQoGgxozxhtnZwfLjlT4n7ipXjJqUs/99RLJDdc98zFmzBhNnTp1yNtlZ2cb22TUo54XatpdzyvPePQaTnZ4/X/itnrJqEk9d9cbbm544yENAABIGTQfAADAKM81H36/X0899ZT8fj/1qOe6esmomYy/0WtS/X/C/Zx6bq93Lte94BQAAKQ2zz3zAQAAvI3mAwAAGEXzAQAAjKL5AAAARtF8AAAAo2g+AACAUTQfAADAKJoPAABg1P8Aj0aXVN9f3wgAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (ax1, ax2) = plt.subplots(1, 2)\n", "ax1.spy(D0.diags(), markersize=2, color='r')\n", "ax2.spy(D1.diags(), markersize=2, color='b')" ] }, { "cell_type": "markdown", "id": "f9c99379", "metadata": { "editable": true }, "source": [ "So $(d^{(0)}_{kj})$ now contains two nonzero diagonals, whereas $(d^{(1)}_{kj})$ is\n", "a matrix with one single diagonal. There is no longer a `full` differentiation\n", "matrix, and we can easily perform this projection for millions of degrees of freedom.\n", "What about $(d^{(2)}_{kj})$? We can now use biharmonic test functions that\n", "satisfy four boundary conditions in the space $B_N = \\{v \\in S_N | v(\\pm 1) = v'(\\pm 1) =0\\}$,\n", "and continue in a similar fashion:" ] }, { "cell_type": "code", "execution_count": 8, "id": "6a5b8db9", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:58.019640Z", "iopub.status.busy": "2024-02-19T13:48:58.019430Z", "iopub.status.idle": "2024-02-19T13:48:58.040659Z", "shell.execute_reply": "2024-02-19T13:48:58.040139Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-6.16153958e-12 -2.01002713e+00 -5.88073092e-12 1.53351702e-01\n", " -2.35856677e-12 -6.03926282e-03 -8.05991953e-13 1.47222318e-04\n", " -2.01978482e-13 -2.61304093e-06 -7.57306469e-29 2.71050543e-20\n", " -1.26217745e-29 2.11758237e-22 0.00000000e+00 0.00000000e+00]\n" ] } ], "source": [ "BN = FunctionSpace(N+4, 'C', bc=(0, 0, 0, 0))\n", "v = TestFunction(BN)\n", "D0 = inner(u, v)\n", "D2 = inner(Dx(u, 0, 2), v)\n", "uN22 = Function(SN)\n", "uN22 = D0.solve(D2.matvec(uN, w0), uN22)\n", "print(uN22-uN2)" ] }, { "cell_type": "markdown", "id": "83944e9d", "metadata": { "editable": true }, "source": [ "We get that `uN22 = uN2`, so the Chebyshev-Petrov-Galerkin projection works. The matrices involved are now" ] }, { "cell_type": "code", "execution_count": 9, "id": "138b71ea", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:58.043558Z", "iopub.status.busy": "2024-02-19T13:48:58.043348Z", "iopub.status.idle": "2024-02-19T13:48:58.349604Z", "shell.execute_reply": "2024-02-19T13:48:58.348943Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAESCAYAAABdK7eSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfPElEQVR4nO3df2zU9R3H8dfR2uNH2nMtgXrjCjUhokURqRoFFaJ2QYLiok7RyvQfyYpQaxh06mRu9EQ3hlsHpP6hLgTlj0llZkY7h0WDP0pL1bkNZHa0ik3HYu4Krge03/3R9EZpKb32+/3cfb/3fCSX5r793r0/Te/9zvvue9/v22dZliUAAABDxiR7AQAAIL3QfAAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMIrmAwAAGEXzAQAAjKL5AAAARtF8AAAAo1zVfGzevFmFhYUaO3as5syZo3fffdeROOFwWFdeeaWys7M1adIkLVmyRAcOHHAk1tni+3w+lZeXOxrnq6++0n333ae8vDyNHz9el19+uRobGx2JderUKT3++OMqLCzUuHHjdOGFF+qpp55ST0+PLc+/Z88eLV68WMFgUD6fT7W1tf1+b1mW1q1bp2AwqHHjxmn+/Pn67LPPHIl38uRJrVmzRpdeeqkmTJigYDCo+++/X0eOHHEk3pkeeugh+Xw+bdq0acTxvIbaYR8v1Q2J2nE6k7XDNc3Hjh07VF5erscee0z79+/Xddddp4ULF6q1tdX2WPX19SorK9MHH3yguro6nTp1SiUlJTp+/Ljtsc7U0NCgmpoaXXbZZY7G+eabbzR37lydd955euONN/S3v/1Nv/rVr3T++ec7Em/Dhg3aunWrqqur9fe//13PPPOMnn32Wf32t7+15fmPHz+uWbNmqbq6etDfP/PMM9q4caOqq6vV0NCg/Px83Xzzzers7LQ93rfffqumpiY98cQTampq0quvvqqDBw/q1ltvHVGsc8U7XW1trT788EMFg8ERx/Iaaod9vFY3JGpHH+O1w3KJq666ylq+fHm/bTNmzLDWrl3reOyOjg5LklVfX+9onM7OTmv69OlWXV2ddcMNN1irVq1yLNaaNWusefPmOfb8Z1q0aJH14IMP9tv2/e9/37rvvvtsjyXJ2rlzZ/x+T0+PlZ+fbz399NPxbV1dXVYgELC2bt1qe7zBfPTRR5Yk6/Dhw47F+/LLL63vfve71l//+ldr6tSp1q9//etRx/ICaod9vFw3LIvaYbJ2uOKTjxMnTqixsVElJSX9tpeUlGjv3r2Ox49EIpKk3NxcR+OUlZVp0aJFuummmxyNI0m7du1ScXGx7rzzTk2aNEmzZ8/W888/71i8efPm6e2339bBgwclSR9//LHee+893XLLLY7F7NPS0qL29vZ+rx+/368bbrjByOtH6n0N+Xw+x94h9vT0qLS0VKtXr1ZRUZEjMdyI2mGvdKobErXDSZnGIo3C0aNH1d3drcmTJ/fbPnnyZLW3tzsa27IsVVRUaN68eZo5c6ZjcV555RU1NTWpoaHBsRin++KLL7RlyxZVVFToJz/5iT766COtXLlSfr9f999/v+3x1qxZo0gkohkzZigjI0Pd3d1av3697rnnHttjnanvNTLY6+fw4cOOx+/q6tLatWu1dOlS5eTkOBJjw4YNyszM1MqVKx15freidtgrneqGRO1wkiuajz4+n6/ffcuyBmyz24oVK/TJJ5/ovffecyxGW1ubVq1apbfeektjx451LM7penp6VFxcrKqqKknS7Nmz9dlnn2nLli2OFJEdO3Zo27Zt2r59u4qKitTc3Kzy8nIFg0EtW7bM9niDScbr5+TJk7r77rvV09OjzZs3OxKjsbFRzz33nJqamhz/e9yK2mGPdKwbErXDCa447DJx4kRlZGQMeKfS0dExoCO108MPP6xdu3Zp9+7dmjJlimNxGhsb1dHRoTlz5igzM1OZmZmqr6/Xb37zG2VmZqq7u9v2mBdccIEuueSSftsuvvhiR76EJ0mrV6/W2rVrdffdd+vSSy9VaWmpHnnkEYXDYUfinS4/P1+SjL9+Tp48qbvuukstLS2qq6tz7J3Lu+++q46ODhUUFMRfP4cPH9ajjz6qadOmORLTLagd9taOdKobErXDSa5oPrKysjRnzhzV1dX1215XV6drr73W9niWZWnFihV69dVX9Ze//EWFhYW2xzjdjTfeqE8//VTNzc3xW3Fxse699141NzcrIyPD9phz584dcArgwYMHNXXqVNtjSb3f4h4zpv/LLSMjw9ZT5s6msLBQ+fn5/V4/J06cUH19vSOvH+n/xePzzz/Xn//8Z+Xl5TkSR5JKS0v1ySef9Hv9BINBrV69Wm+++aZjcd2A2mFv7UinuiFRO5zkmsMuFRUVKi0tVXFxsa655hrV1NSotbVVy5cvtz1WWVmZtm/frtdee03Z2dnxrjcQCGjcuHG2x8vOzh5wTHjChAnKy8tz7FjxI488omuvvVZVVVW666679NFHH6mmpkY1NTWOxFu8eLHWr1+vgoICFRUVaf/+/dq4caMefPBBW57/2LFjOnToUPx+S0uLmpublZubq4KCApWXl6uqqkrTp0/X9OnTVVVVpfHjx2vp0qW2xwsGg7rjjjvU1NSk119/Xd3d3fHXUG5urrKysmz/+84sUOedd57y8/N10UUXjejv8xJqh328VjckakfSaofj59PY6He/+501depUKysry7riiiscO31N0qC3F154wZF4g3H6VFvLsqw//vGP1syZMy2/32/NmDHDqqmpcSxWNBq1Vq1aZRUUFFhjx461LrzwQuuxxx6zYrGYLc+/e/fuQf9ny5Ytsyyr95S5J5980srPz7f8fr91/fXXW59++qkj8VpaWs76Gtq9e7cjf9+ZONW2P2qHfbxUNyyL2nEmU7XDZ1mWZWczAwAAMBRXfOcDAAB4B80HAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRrms+YrGY1q1bp1gsRjzipVy8ZMRMxt/oNl7/n/A6J16qxzuT6y4yFo1GFQgEFIlEHBu2QzziuSlmMv5Gt/H6/4TXOfFSPd6ZXPfJBwAAcDeaDwAAYFTKTbXt6enRkSNHlJ2dLZ/PN+D30Wi030+nEY94qR7TiXiWZamzs1PBYHDASPNUNVTt8ML/JJXiJSMm8VI/XiJ1I+W+8/Hll18qFAolexkAJLW1tWnKlCnJXsawUDuA1DCcupFyn3xkZ2dL6l08X54DkiMajSoUCsXz0Q2oHUByJVI3HGs+Nm/erGeffVZff/21ioqKtGnTJl133XXnfFzfx6U5OTkUECDJBjv06aSR1g2J2gGkiuHUDUcO5u7YsUPl5eV67LHHtH//fl133XVauHChWltbnQgHwAOoG0D6cOQ7H1dffbWuuOIKbdmyJb7t4osv1pIlSxQOh4d8bLLPPQaQnDwcTd2QqB1AsiWSg7Z/8nHixAk1NjaqpKSk3/aSkhLt3bt39AFaW6Wmpt6fJng9HpACHK8bhpHGwNBs/87H0aNH1d3drcmTJ/fbPnnyZLW3tw/YPxaL9bu865Cn/bS2ShddJHV1SWPHSgcOSAUFtq097eIBKSLRuiElWDsMIo2Bc3PsBP4zv3BiWdagX0IJh8MKBALx25Cnyh092pvRUu/Po0ftXHL6xQNSzHDrhpRg7TCINAbOzfbmY+LEicrIyBjwbqWjo2PAuxpJqqysVCQSid/a2tqGevLetxJS78+JE+1cevrFA1JEonVDSrB2GEQaA+dm+2GXrKwszZkzR3V1dbr99tvj2+vq6nTbbbcN2N/v98vv9w/vyQsKej/DPHq0N6Od/izT6/GAFJFo3ZASrB0GkcbAuTlynY+KigqVlpaquLhY11xzjWpqatTa2qrly5eP/skLCsxms9fjASnC0bphGGkMDM2R5uMHP/iB/vOf/+ipp57S119/rZkzZ+pPf/qTpk6d6kQ4AB5A3QDSR8rNduFcfSD53JiHblwz4CVJvc4HAADAUGg+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAY5ch1PlJSa6vZSw56PR4A25HGSBfp0Xx4fTotYzQB1yONkU7S47CL16fTMkYTcD3SGOkkPZoPr0+nZYwm4HqkMdJJehx28fp0WsZoAq5HGiOdpEfzIXl/Oi1jNAHXI42RLtLjsAsAAEgZNB8AAMAomg8AAGAUzQcAADCK5gMAABhle/MRDod15ZVXKjs7W5MmTdKSJUt04MABu8MA8BDqBpBebG8+6uvrVVZWpg8++EB1dXU6deqUSkpKdPz4cbtDAfAI6gaQXnyWZVlOBvj3v/+tSZMmqb6+Xtdff/05949GowoEAopEIsrJyXFyaQDOItl5mGjdkJK/ZiDdJZKDjl9kLBKJSJJyc3OdDuUMr0+nZYwmUpDr64ZhpDHcxtHmw7IsVVRUaN68eZo5c+ag+8RiMcVisfj9aDTq5JIS4/XptIzRRAoaTt2QUrx2GEQaw40cPdtlxYoV+uSTT/Tyyy+fdZ9wOKxAIBC/hUIhJ5eUGK9Pp2WMJlLQcOqGlOK1wyDSGG7kWPPx8MMPa9euXdq9e7emTJly1v0qKysViUTit7a2NqeWlDivT6dljCZSzHDrhpTitcMg0hhuZPthF8uy9PDDD2vnzp165513VFhYOOT+fr9ffr/f7mXYw+vTaRmjiRSRaN2QUrx2GEQaw41sbz7Kysq0fft2vfbaa8rOzlZ7e7skKRAIaNy4cXaHc57Xp9MyRhMpwHN1wzDSGG5j+6m2Pp9v0O0vvPCCfvjDH57z8ZwuBySf6Twcbd2QqB1AsiX1VFuHLxsCwIOoG0B6YbYLAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRjk+1TVten07LGE3AE0hlJAPNhxO8Pp2WMZqAJ5DKSBYOuzjB69NpGaMJeAKpjGSh+XCC16fTMkYT8ARSGcnCYRcneH06LWM0AU8glZEsNB9O8fp0WsZoAp5AKiMZOOwCAACMovkAAABG0XwAAACjaD4AAIBRNB8AAMAomg8AAGCU481HOByWz+dTeXm506EAeAR1A/A2R5uPhoYG1dTU6LLLLnMyDAAPoW4A3udY83Hs2DHde++9ev755/Wd73zHqTDe09oqNTX1/iQe0gx1wz1IZYyGY81HWVmZFi1apJtuumnI/WKxmKLRaL9b2uobMTlnTu9Pp7Pa6/HgOsOtGxK1I5lIZYyWI83HK6+8oqamJoXD4XPuGw6HFQgE4rdQKOTEktzB69NpGaGJISRSNyRqRzKRyhgt25uPtrY2rVq1Stu2bdPYvnGJQ6isrFQkEonf2tra7F6Se3h9Oi0jNHEWidYNidqRTKQyRstnWZZl5xPW1tbq9ttvV0ZGRnxbd3e3fD6fxowZo1gs1u93Z4pGowoEAopEIsrJybFzae7Q2mp2xKTX42FETOfhaOuGRO0wjVTGmRLJQdun2t5444369NNP+2174IEHNGPGDK1Zs+acBSTteX06LSM0MQjqhvuQyhgN25uP7OxszZw5s9+2CRMmKC8vb8B2AJCoG0C64QqnAADAKNs/+RjMO++8YyIMAA+hbgDexScfAADAKJoPAABgFM0HAAAwiuYDAAAYRfMBAACMMnK2Cwzw+pVKuZwi4AmkMiSaD2/oGzHZ1dU7aOHAAWez2uvxADiCVEYfDrt4gden0zJCE/AEUhl9aD68wOvTaRmhCXgCqYw+HHbxgoKC3s8vTR1I9Xo8AI4gldGH5sMrvD6dlhGagCeQypA47AIAAAyj+QAAAEbRfAAAAKNoPgAAgFE0HwAAwCiaDwAAYJQjzcdXX32l++67T3l5eRo/frwuv/xyNTY2OhEKgEdQN4D0Yft1Pr755hvNnTtXCxYs0BtvvKFJkybpn//8p84//3y7QwHwCOoGkF5sbz42bNigUCikF154Ib5t2rRpdoeBXbw+nZYRmq5A3cC5kMreYvthl127dqm4uFh33nmnJk2apNmzZ+v5558/6/6xWEzRaLTfDYb0jZicM6f3Z2sr8ZAUidYNidqRTkhl77G9+fjiiy+0ZcsWTZ8+XW+++aaWL1+ulStX6ve///2g+4fDYQUCgfgtFArZvSScjden0zJC0zUSrRsStSOdkMre47Msy7LzCbOyslRcXKy9e/fGt61cuVINDQ16//33B+wfi8UUi8Xi96PRqEKhkCKRiHJycuxcGs7U93aiq6t3xOSBA85+nun1eB4SjUYVCASM5WGidUOidqQTUtkdEqkbtn/n44ILLtAll1zSb9vFF1+sP/zhD4Pu7/f75ff77V4GhsPr02kZoekaidYNidqRTkhl77G9+Zg7d64OHDjQb9vBgwc1depUu0PBDl6fTssITVegbuBcSGVvsf07H4888og++OADVVVV6dChQ9q+fbtqampUVlZmdygAHkHdANKL7c3HlVdeqZ07d+rll1/WzJkz9fOf/1ybNm3Svffea3coAB5B3QDSi+1fOB0t0190AzCQG/PQjWsGvCSRHGS2CwAAMIrmAwAAGEXzAQAAjKL5AAAARtF8AAAAo2y/yBjSRDpMp2WMJuB6pHFqovlA4tJhRgvDJADXI41TF4ddkLh0mE7LGE3A9Ujj1EXzgcRNnNj7NkLq/TlxorfiJSsmAFuRxqmLwy5IXDpMp2WMJuB6pHHqovnAyKTDdFrGaAKuRxqnJg67AAAAo2g+AACAUTQfAADAKJoPAABgFM0HAAAwiuYDAAAYZXvzcerUKT3++OMqLCzUuHHjdOGFF+qpp55ST0+P3aEAeAR1A0gvtl/nY8OGDdq6dateeuklFRUVad++fXrggQcUCAS0atUqu8MB8ADqBpBebG8+3n//fd12221atGiRJGnatGl6+eWXtW/fPrtDwY2YhotBUDeQakhjZ9l+2GXevHl6++23dfDgQUnSxx9/rPfee0+33HLLoPvHYjFFo9F+N3hU34jJOXN6f7a2eitesmJ6QKJ1Q6J2wDmksfNs/+RjzZo1ikQimjFjhjIyMtTd3a3169frnnvuGXT/cDisn/3sZ3YvA6losBGTTr6lMB0vWTE9ING6IVE74BzS2Hm2f/KxY8cObdu2Tdu3b1dTU5Neeukl/fKXv9RLL7006P6VlZWKRCLxW1tbm91LQqpgGi7OItG6IVE74BzS2Hk+y7IsO58wFApp7dq1Kisri2/7xS9+oW3btukf//jHOR8fjUYVCAQUiUSUk5Nj59KQCvjOhyuYzsPR1g2J2gF7eSCNjUskB20/7PLtt99qzJj+H6hkZGRwyhx6MQ0Xg6BuINWQxs6yvflYvHix1q9fr4KCAhUVFWn//v3auHGjHnzwQbtDAfAI6gaQXmw/7NLZ2aknnnhCO3fuVEdHh4LBoO655x799Kc/VVZW1jkfz0enQPKZzsPR1g2J2gEkWyI5aHvzMVoUECD53JiHblwz4CWJ5CCzXQAAgFE0HwAAwCiaDwAAYBTNBwAAMIrmAwAAGEXzAQAAjLL9ImOAI9LhMulczxlwPdJ4eGg+kPr65lt3dfVOeTpwwPmsNh0zGX8jAFuRxsPHYRekvsHmW3stZjL+RgC2Io2Hj+YDqS8Z861Nx2SGN+B6pPHwcdgFqa+goPfzS5MHUk3HTMbfCMBWpPHw0XzAHZIx39p0TGZ4A65HGg8Ph10AAIBRNB8AAMAomg8AAGAUzQcAADCK5gMAABiVcPOxZ88eLV68WMFgUD6fT7W1tf1+b1mW1q1bp2AwqHHjxmn+/Pn67LPP7FovABeibgA4XcLNx/HjxzVr1ixVV1cP+vtnnnlGGzduVHV1tRoaGpSfn6+bb75ZnZ2do14sAHeibgDoxxoFSdbOnTvj93t6eqz8/Hzr6aefjm/r6uqyAoGAtXXr1mE9ZyQSsSRZkUhkNEsDMApO5qETdcOyqB1AsiWSg7Z+56OlpUXt7e0qKSmJb/P7/brhhhu0d+9eO0MBw9PaKjU19f70asxk/I02om4Ark/jhNl6hdP29nZJ0uTJk/ttnzx5sg4fPjzoY2KxmGKxWPx+NBq1c0lIZ0zDdYWR1A2J2gHv8EAaJ8yRs118Pl+/+5ZlDdjWJxwOKxAIxG+hUMiJJSEdMQ3XVRKpGxK1A97hoTQeNlubj/z8fEn/fyfTp6OjY8C7mj6VlZWKRCLxW1tbm51LQjpjGq4rjKRuSNQOeIcH0jhhth52KSwsVH5+vurq6jR79mxJ0okTJ1RfX68NGzYM+hi/3y+/32/nMoBeTMN1hZHUDYnaAe/wQBonLOHm49ixYzp06FD8fktLi5qbm5Wbm6uCggKVl5erqqpK06dP1/Tp01VVVaXx48dr6dKlti4cGBam4aYE6gYwNBeksa0Sbj727dunBQsWxO9XVFRIkpYtW6YXX3xRP/7xj/Xf//5XP/rRj/TNN9/o6quv1ltvvaXs7Gz7Vg3AVagbAE7nsyzLSvYiTheNRhUIBBSJRJSTk5Ps5QBpyY156MY1A16SSA4y2wUAABhF8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFE0HwAAwChbr3AKeEprq9lLDno9HgDbuTWNaT6AwXh9Om06jtEEPMbNacxhF2AwXp9Om45jNAGPcXMa03wAg/H6dNp0HKMJeIyb05jDLsBgvD6dNh3HaAIe4+Y0pvkAzsbr02nTbYwm4EFuTWMOuwAAAKNoPgAAgFE0HwAAwCiaDwAAYBTNBwAAMCrh5mPPnj1avHixgsGgfD6famtr4787efKk1qxZo0svvVQTJkxQMBjU/fffryNHjti5ZgAuQ90AcLqEm4/jx49r1qxZqq6uHvC7b7/9Vk1NTXriiSfU1NSkV199VQcPHtStt95qy2IBuBN1A8DpfJZlWSN+sM+nnTt3asmSJWfdp6GhQVdddZUOHz6sgmGcjByNRhUIBBSJRJSTkzPSpQEYBSfz0Im6IVE7gGRLJAcdv8hYJBKRz+fT+eef73QoIDUwnXbUqBuAs5JdNhxtPrq6urR27VotXbr0rF1QLBZTLBaL349Go04uCXAW02lHbTh1Q6J2ACOVCmXDsbNdTp48qbvvvls9PT3avHnzWfcLh8MKBALxWygUcmpJgPOYTjsqw60bErUDGKlUKBuONB8nT57UXXfdpZaWFtXV1Q357qWyslKRSCR+a2trc2JJgBlMpx2xROqGRO0ARioVyobth136Csjnn3+u3bt3Ky8vb8j9/X6//H6/3csAkoPptCOSaN2QqB3ASKVC2Ui4+Th27JgOHToUv9/S0qLm5mbl5uYqGAzqjjvuUFNTk15//XV1d3ervb1dkpSbm6usrCz7Vg6kKqbTDkDdAFJLsstGwqfavvPOO1qwYMGA7cuWLdO6detUWFg46ON2796t+fPnn/P5OV0OSD6789DpuiFRO4Bkc/RU2/nz52uofmUUlw0B4FHUDQCnY7YLAAAwiuYDAAAYRfMBAACMovkAAABG0XwAAACjaD4AAIBRjk+1BTBMTMMFkOLsKhs0H0AqYBougBRnZ9ngsAuQCpiGCyDF2Vk2aD6AVMA0XAApzs6ywWEXIBUwDRdAirOzbNB8AKmCabgAUpxdZYPDLgAAwCiaDwAAYBTNBwAAMIrmAwAAGEXzAQAAjKL5AAAARiXcfOzZs0eLFy9WMBiUz+dTbW3tWfd96KGH5PP5tGnTplEsEYDbUTcAnC7h5uP48eOaNWuWqqurh9yvtrZWH374oYLB4IgXB8AbqBsATpfwRcYWLlyohQsXDrnPV199pRUrVujNN9/UokWLRrw4AMPggum01A0Ap7P9Cqc9PT0qLS3V6tWrVVRUdM79Y7GYYrFY/H40GrV7SYB3eWQ6baJ1Q6J2AG5m+xdON2zYoMzMTK1cuXJY+4fDYQUCgfgtFArZvSTAuzwynTbRuiFROwA3s7X5aGxs1HPPPacXX3xRPp9vWI+prKxUJBKJ39ra2uxcEuBtHphOO5K6IVE7ADeztfl499131dHRoYKCAmVmZiozM1OHDx/Wo48+qmnTpg36GL/fr5ycnH43AMPUN2aysdG1h1xGUjckagfgZrZ+56O0tFQ33XRTv23f+973VFpaqgceeMDOUAD6uHw6LXUDSD8JNx/Hjh3ToUOH4vdbWlrU3Nys3NxcFRQUKC8vr9/+5513nvLz83XRRReNfrUAXIm6AeB0CTcf+/bt04IFC+L3KyoqJEnLli3Tiy++aNvCAHgHdQPA6RJuPubPny/Lsoa9/7/+9a9EQwDwGOoGgNMx2wUAABhF8wEAAIyi+QAAAEbRfAAAAKNoPgAAgFG2D5YD4BIumIYLwJtoPoB05JFpuADcicMuQDryyDRcAO5E8wGkIw9MwwXgXhx2AdJR3zRcvvMBIAloPoB05fJpuADci8MuAADAKJoPAABgVModdumbfBmNRpO8EiB99eVfIpNok43aASRXInUj5ZqPzs5OSVIoFErySgB0dnYqEAgkexnDQu0AUsNw6obPSrG3Nj09PTpy5Iiys7Pl8/kG/D4ajSoUCqmtrU05OTmOr4d4xEv1mE7EsyxLnZ2dCgaDGjPGHUdnh6odXvifpFK8ZMQkXurHS6RupNwnH2PGjNGUKVPOuV9OTo6xJCMe8dwQ0+54bvnEo89waofb/yepFi8ZMYmX2vGGWzfc8ZYGAAB4Bs0HAAAwynXNh9/v15NPPim/30884qVcvGTETMbf6DZe/5/wOideqsc7U8p94RQAAHib6z75AAAA7kbzAQAAjKL5AAAARtF8AAAAo2g+AACAUTQfAADAKJoPAABgFM0HAAAw6n9FUopHGPWkLgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (ax1, ax2) = plt.subplots(1, 2)\n", "ax1.spy(D0.diags(), markersize=2, color='r')\n", "ax2.spy(D2.diags(), markersize=2, color='b')" ] }, { "cell_type": "markdown", "id": "99a670bc", "metadata": { "editable": true }, "source": [ "So there are now three nonzero diagonals in $(d^{(0)}_{kj})$, whereas the differentiation matrix\n", "$(d^{(2)}_{kj})$ contains only one nonzero diagonal.\n", "\n", "Why does this work so well? The Chebyshev polynomials and their derivatives satisfy the following orthogonality relation" ] }, { "cell_type": "markdown", "id": "e9c66875", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eq:orthon} \\tag{6}\n", " \\int_{-1}^{1} T^{(n)}_j T^{(n)}_k \\omega^{n-1/2} dx = \\alpha^{n}_k \\delta_{kj}, \\quad \\text{for}\\, n \\ge 0,\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ab4aeed8", "metadata": { "editable": true }, "source": [ "where" ] }, { "cell_type": "markdown", "id": "2e7f11e2", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " \\alpha^n_k = \\frac{c_{k+n}\\pi k (k+n-1)!}{2(k-n)!}.\n", "\\label{_auto3} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "ec8e259a", "metadata": { "editable": true }, "source": [ "So when we choose a test function $\\omega^n T^{(n)}_k$ and a trial function $T_j$,\n", "we get the diagonal differentiation matrix" ] }, { "cell_type": "markdown", "id": "0a9845f8", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", " d^{(n)}_{kj} = \\int_{-1}^{1} T^{(n)}_j (\\omega^n T^{(n)}_k) \\omega^{-1/2} dx = \\alpha^{n}_k \\delta_{kj}, \\quad \\text{for}\\, n \\ge 0.\n", "\\label{_auto4} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1faaf249", "metadata": { "editable": true }, "source": [ "The two chosen test functions above are both proportional to $\\omega^n T^{(n)}_k$. More precisely,\n", "$T_k-T_{k+2} = \\frac{2}{k+1} \\omega T^{(1)}_{k+1}$ and the biharmonic test function is\n", "$T_k-\\frac{2(k+2)}{k+3}T_{k+2} + \\frac{k+1}{k+3}T_{k+4} = \\frac{4 \\omega^2T^{(2)}_{k+2}}{(k+2)(k+3)}$.\n", "Using these very specific test functions correponds closely to using the Chebyshev\n", "recursion formulas that are found in most textbooks. Here they are adapted to\n", "a Chebyshev-Petrov-Galerkin method, where we simply choose test and trial functions and everything\n", "else falls into place in a few lines of code." ] }, { "cell_type": "markdown", "id": "505a7c63", "metadata": { "editable": true }, "source": [ "## Recursion\n", "\n", "Let us for completion show how to\n", "find $\\hat{u}_N^{(1)}$ with a recursive approach. The Chebyshev polynomials\n", "satisfy" ] }, { "cell_type": "markdown", "id": "e01fad08", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "2T_k = \\frac{1}{k+1}T'_{k+1}- \\frac{1}{k-1} T'_{k-1}, \\quad k \\ge 1.\n", "\\label{eq:Trec1} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "a5e54439", "metadata": { "editable": true }, "source": [ "By using this and setting $u' = \\sum_{k=0}^{\\infty} \\hat{u}^{(1)}_k T_k = \\sum_{k=0}^{\\infty} \\hat{u}_k T'_k$\n", "we get" ] }, { "cell_type": "markdown", "id": "c6be810e", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "2k\\hat{u}_k = c_{k-1}\\hat{u}^{(1)}_{k-1} - \\hat{u}^{(1)}_{k+1}, \\quad k \\ge 1.\n", "\\label{eq:Trec2} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "235acf35", "metadata": { "editable": true }, "source": [ "Using this recursion on a discrete series, together with $\\hat{u}^{(1)}_{N} = \\hat{u}^{(1)}_{N-1} = 0$, we get\n", "(see [[canuto]](#canuto) Eq. (2.4.24))" ] }, { "cell_type": "markdown", "id": "a3b2391d", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "c_k \\hat{u}^{(1)}_k = \\hat{u}^{(1)}_{k+2} + 2(k+1)\\hat{u}_{k+1}, \\quad 0 \\le k \\le N-2,\n", "\\label{eq:Trec3} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "acbd18e5", "metadata": { "editable": true }, "source": [ "which is easily implemented in a (slow) for-loop" ] }, { "cell_type": "code", "execution_count": 10, "id": "b62bbb21", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:58.352963Z", "iopub.status.busy": "2024-02-19T13:48:58.352760Z", "iopub.status.idle": "2024-02-19T13:48:58.357355Z", "shell.execute_reply": "2024-02-19T13:48:58.356774Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.00000000e+00 0.00000000e+00 0.00000000e+00 7.88860905e-31\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.88860905e-31\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n" ] } ], "source": [ "f1 = np.zeros(N+1)\n", "ck = np.ones(N); ck[0] = 2\n", "for k in range(N-2, -1, -1):\n", " f1[k] = (f1[k+2]+2*(k+1)*uN[k+1])/ck[k]\n", "print(f1[:-1]-uN1)" ] }, { "cell_type": "markdown", "id": "4dcda2f1", "metadata": { "editable": true }, "source": [ "which evidently is exactly the same result. It turns out that this is not strange. If we multiply\n", "([11](#eq:Trec3)) by $\\pi/2$, rearrange a little bit and move to a matrix form we get" ] }, { "cell_type": "markdown", "id": "eaeaa2a8", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "c_k \\pi/2 \\hat{u}^{(1)}_k - \\pi/2 \\hat{u}^{(1)}_{k+2} = (k+1)\\pi \\hat{u}_{k+1} \n", "\\label{_auto5} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "1eddc342", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\underbrace{(c_k \\pi/2 \\delta_{kj} - \\pi/2 \\delta_{k, j-2})}_{(D^0)_{kj}} \\hat{u}^{(1)}_j = \\underbrace{(k+1)\\pi \\delta_{k,j-1}}_{(D^1)_{kj}} \\hat{u}_{j} \n", "\\label{_auto6} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "f3e68c7c", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "D^0 \\boldsymbol{\\hat{u}} = D^1 \\boldsymbol{\\hat{u}} \n", "\\label{_auto7} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "e4d34134", "metadata": { "editable": true }, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\boldsymbol{\\hat{u}^{(1)}} = (D^0)^{-1} D^1 \\boldsymbol{\\hat{u}}\n", "\\label{_auto8} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "b2e1a8af", "metadata": { "editable": true }, "source": [ "which is exactly how $\\boldsymbol{\\hat{u}^{(1)}}$ was computed above with the Chebyshev-Petrov-Galerkin approach\n", "(compare with the code line `uN11 = D0.solve(D1.matvec(uN, w0), uN11)`). Not convinced? Check that the matrices\n", "`D0` and `D1` are truly as stated above. The matrices below are printed as dictionaries with diagonal\n", "number as key (main is 0, first upper is 1 etc) and diagonal values as values:" ] }, { "cell_type": "code", "execution_count": 11, "id": "f76f13c5", "metadata": { "collapsed": false, "editable": true, "execution": { "iopub.execute_input": "2024-02-19T13:48:58.360072Z", "iopub.status.busy": "2024-02-19T13:48:58.359900Z", "iopub.status.idle": "2024-02-19T13:48:58.366474Z", "shell.execute_reply": "2024-02-19T13:48:58.365608Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{0: array([3.14159265, 1.57079633, 1.57079633, 1.57079633, 1.57079633,\n", " 1.57079633, 1.57079633, 1.57079633, 1.57079633, 1.57079633,\n", " 1.57079633, 1.57079633, 1.57079633, 1.57079633, 1.57079633,\n", " 1.57079633]),\n", " 2: array([-1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633,\n", " -1.57079633, -1.57079633, -1.57079633, -1.57079633, -1.57079633,\n", " -1.57079633, -1.57079633, -1.57079633, -1.57079633])}\n", "{1: array([ 3.14159265, 6.28318531, 9.42477796, 12.56637061, 15.70796327,\n", " 18.84955592, 21.99114858, 25.13274123, 28.27433388, 31.41592654,\n", " 34.55751919, 37.69911184, 40.8407045 , 43.98229715, 47.1238898 ,\n", " 50.26548246])}\n" ] } ], "source": [ "import pprint\n", "DN = FunctionSpace(N+2, 'C', bc=(0, 0))\n", "v = TestFunction(DN)\n", "D0 = inner(u, v)\n", "D1 = inner(Dx(u, 0, 1), v)\n", "pprint.pprint(dict(D0))\n", "pprint.pprint(dict(D1))" ] }, { "cell_type": "markdown", "id": "a43c5ad2", "metadata": { "editable": true }, "source": [ "In conclusion, we have shown that we can use an efficient Chebyshev-Petrov-Galerkin approach to obtain\n", "the discrete Chebyshev coefficients for the derivatives\n", "of a function. By inspection, it turns out that this approach is identical to the common methods based on\n", "well-known Chebyshev recursion formulas.\n", "\n", "\n", "\n", "1.
**C. Canuto, M. Hussaini, A. Quarteroni and J. T. A.**. *Spectral Methods in Fluid Dynamics*, *Scientific Computation*, Springer, 2012." ] } ], "metadata": { "kernelspec": { "display_name": "shenfun24", "language": "python3", "name": "shenfun24" }, "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.0" } }, "nbformat": 4, "nbformat_minor": 5 }