{ "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": [ "