{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# TD2: Spatial Discretization\n", "\n", "With spatial discretization, the goal is to compute all values $f(x_0), ..., f(x_N)$ of a discretized function $f$ simultaneously. Spatial discretization is therefore opposed to temporal discretization for which the computation of $f(t_{i+1})$ relies on the previous values $f(t_0),..., f(t_i)$.\n", "\n", "In this TD we will not be interested in solving a differential equation with spatial discretization as we first need to get familiar with methods to approximate a function $f$, and we will focus here on **spectral methods**. We recall that spectral methods rely on the expansion of a function $f$ on a basis of functions $\\{\\phi_k\\}_k$ (in a function space we do not specify at this point). The coefficients of the expansion will be noted $\\hat{f}_k$. Thus we have:\n", "$$\n", " f(x) = \\sum_{k}\\hat{f}_k\\phi_k(x)\n", "$$\n", "\n", "The choice of the basis $\\{\\phi_k\\}_k$ is obviously crucial for the efficiency of a numerical calculation. We will make sure for a practical case that the chosen basis meets the following requirements:\n", "* **Convergence:** the series should converge rapidly enough to the function $f$\n", "* **Transformation:** the computation of the coefficients $\\hat{f}_k$ should be numerically fast, as well as the inverse transformation\n", "* **Differentiation:** it should be easy to compute the derivative of $f$ through its expansion, meaning we can determine $\\widehat{f^\\prime}_k$ easily\n", "\n", "As usual, we will need the following packages:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: Fourier and Chebyshev interpolation\n", "\n", "In spectral methods, given a set of points $\\{x_k\\}_k$ we define the interpolant $\\mathbb{I}_N[f]$ of a function $f$ based on the series representation of $f$ on the basis $\\{\\phi_k\\}_k$. Thus, we approximate $f$ by a truncated series:\n", "$$\n", " \\mathbb{I}_N[f] = \\sum_{k = 0}^{N}\\hat{f}_k\\phi_k\n", "$$\n", "(depending on the chosen basis, we can have slightly different definitions of $\\mathbb{I}_N[f]$). This function in an interpolation of $f$, and should therefore satisfy $\\forall k, f(x_k) = \\mathbb{I}_N[f](x_k)$.\n", "\n", "### Fourier basis\n", "\n", "First let us consider a Fourier basis $\\phi_k(x) = e^{i k x}$ in the function space $L^2([0,2\\pi])$. We choose the spatial domain $[0, 2\\pi]$ for convenience in the expression of the following formula, but any scaling would make the following transposable to any compact domain $[a,b]$. We also define a set of $N$ regularly spaced points in this domain $\\{x_m = \\frac{2\\pi m}{N}, m = 0, ..., N - 1\\}$. Now for a function $f$, we define its approximation by the interpolant:\n", "$$\n", " \\mathbb{I}_N[f](x) = \\sum_{k = \\left \\lfloor{-N/2}\\right \\rfloor + 1}^{\\left \\lfloor{N/2}\\right \\rfloor}{\\hat{f}_k e^{ikx}}\n", "$$\n", "The determination of the interpolant is directly related to the FFT such as defined in the [*numpy* implementation](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html): the coefficients $\\hat{f}_k$ are those computed by *np.fft.fft (fVals) / N* where *fVals* is the vector of values of the set $\\{f(x_m)\\}_m$ and *N* the length of *fVals*. We will pay attention to the order of the coefficients returned by *np.fft.fft* function. We also recall that the FFT algorithm has a $\\mathcal{O}(N\\log(N))$ complexity.\n", "\n", "**1)** Write a function which computes from a given set of values $\\{f(x_m)\\}_m$ the interpolant coefficients $\\{\\hat{f}_k\\}_k$. Complexity should be a $\\mathcal{O}(N\\log(N))$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**2)** Write a function which computes from an arbitrary set of points and for a given set of coefficients $\\{\\hat{f}_k\\}_k$ the values of $\\mathbb{I}_N[f]$ evaluated at the input points. We will not care about the complexity of the algorithm." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**3)** Compare by plotting on $[0,2\\pi]$ in the following cases the exact function with its interpolant for several values of $N$:\n", "* $f(x) = \\cos(x) + \\sin(2x)$\n", "* $f(x) = \\begin{cases}\n", " \\frac{\\pi}{4} & \\text{if } 0 \\leq x < \\pi \\\\\n", " -\\frac{\\pi}{4} & \\text{if } \\pi \\leq x < 2\\pi\n", " \\end{cases}$\n", "* $f(x) = \\frac1{1+25(x/\\pi - 1)^2}$\n", " \n", "We will pay attention to the emergence of the Gibbs phenomenon." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Chebyshev basis\n", "\n", "To avoid some unwanted numerical behaviors such as Gibbs phenomenon when dealing with discontinuities, we can use another basis of functions: Chebyshev polynomials.\n", "\n", "The Chebyshev polynomial of degree $n$ is the unique polynomial satisfying:\n", "$$\n", " \\forall \\theta \\in \\mathbb{R}, T_n(\\cos(\\theta)) = \\cos(n\\theta)\n", "$$\n", "We can build Chebyshev polynomials with the following recursive definition:\n", "$$\n", " \\begin{cases}\n", " T_0(X) = 1\\\\\n", " T_1(X) = X\\\\\n", " T_{n+2}(X) = 2XT_{n+1}(X) - T_{n}(X) \\text{ with } n \\in \\mathbb{N}\n", " \\end{cases}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**4)** Write a function which computes the $N + 1$ first Chebyshev polynomials for a given $N$, returning the coefficients of the polynomials in a matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**5)** Write a function which computes from an arbitrary set of points and for a given set of coefficients $\\{a_k\\}$ the values of the polynomial $\\sum_k{a_k}X^k$ evaluated at the input points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**6)** Using the previous functions, plot the first Chebyshev polynomials." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a more concrete view on Chebyshev polynomials, we can use them to approximate a function $f$. We will now restrict to a function $f$ which is defined on the spatial domain $[-1,1]$, but as before, this is easily generalizable to a function defined on any compact domain $[a,b]$ through a scaling.\n", "\n", "For a given $N$, we define a grid of the spatial domain as follows:\n", "$$\n", " x_j = \\cos(\\theta_j) \\text{ with } \\theta_j = \\frac{\\pi j}{N} \\text{ and } j=0,...,N\n", "$$\n", "We can note we have one more point on the grid in comparison with the Fourier basis for the same value of N.\n", "We also define the interpolant of $f$ as follows:\n", "$$\n", " \\mathbb{I}_N[f](x) = \\sum_{n=0}^N\\hat{f}_nT_n(x)\n", "$$\n", "\n", "The determination of the interpolant for a given function on the previsouly defined grid also relies on the FFT. More precisely, it can be computed using the [Discrete Cosine Transform](https://en.wikipedia.org/wiki/Discrete_cosine_transform). To save some time, we can use the following function which computes the coefficients of the interpolant from a given set of values $\\{f(x_j)\\}_j$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy as scp\n", "import scipy.fftpack\n", "\n", "def ChebTransform (x):\n", " ret = scp.fftpack.dct (x, type = 1) / (2 * (len (x) - 1))\n", " ret [1:-1] = ret [1:-1] * 2\n", " return ret" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function has a $\\mathcal{O}(N\\log(N))$ complexity." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**7)** Write a function which computes from an arbitrary set of points and for a given set of coefficients $\\{\\hat{f}_n\\}_n$ the values of $\\mathbb{I}_N[f]$ evaluated at the input points. We will not care about the complexity of the algorithm." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**8)** Compare by plotting on $[-1,1]$ in the following cases the exact function with its interpolant for several values of $N$:\n", "* $f(x) = \\cos((x + 1)\\pi) + \\sin(2(x + 1)\\pi)$\n", "* $f(x) = \\begin{cases}\n", " \\frac{\\pi}{4} & \\text{if } -1 \\leq x < 0 \\\\\n", " -\\frac{\\pi}{4} & \\text{if } 0 \\leq x < 1\n", " \\end{cases}$\n", "* $f(x) = 1/(1+25x^2)$\n", "\n", "In which cases Chebyshev interpolation seems to be more suitable in comparison with Fourier interpolation?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: Spectral differentiation\n", "\n", "Spectral representations will be useful to solve a differential equation only if we can easily differentiate them. In this part, we will see how we can do that within a Fourier representation and a Chebyshev representation of a given function $f$.\n", "\n", "### Fourier representation\n", "\n", "Using the same notations as before, we already defined an interpolant $\\mathbb{I}_N[f]$ which approximates $f$. It is easy to differentiate $\\mathbb{I}_N[f]$:\n", "$$\n", " \\mathbb{I}_N[f]^\\prime(x)= \\sum_{k = \\left \\lfloor{-N/2}\\right \\rfloor + 1}^{\\left \\lfloor{N/2}\\right \\rfloor}{ik\\hat{f}_k e^{ikx}}\n", "$$\n", "Now we can wonder if this function is still a good approximation of $f^\\prime$.\n", "\n", "**1)** Write a function which computes the coefficients of $\\mathbb{I}_N[f]^\\prime$ given the set of $\\{\\hat{f}_k\\}_k$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**2)** Using *np.fft.ifft* function ($\\mathcal{O}(N\\log(N))$ complexity), write a function which computes the values $\\{\\mathbb{I}_N[f]^\\prime(x_m)\\}_m$ given the coefficients of $\\mathbb{I}_N[f]^\\prime$. Complexity of the algorithm should be a $\\mathcal{O}(N\\log(N))$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**3)** Compare by plotting the exact value of $f^\\prime$ and the approximation $\\mathbb{I}_N[f]^\\prime$ in the following cases:\n", "* $f(x) = \\cos(x) + \\sin(2x)$\n", "* $f(x) = \\begin{cases}\n", " \\frac{\\pi}{4} & \\text{if } 0 \\leq x < \\pi \\\\\n", " -\\frac{\\pi}{4} & \\text{if } \\pi \\leq x < 2\\pi\n", " \\end{cases}$\n", "* $f(x) = \\frac1{1+25(x/\\pi - 1)^2}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Chebyshev representation\n", "\n", "Exactly on the same model as before, and using the same notations as part 1, we wonder if $\\mathbb{I}_N[f]^\\prime$ is a good approximation of $f$, where:\n", "$$\n", " \\mathbb{I}_N[f]^\\prime(x) = \\sum_{n=0}^N\\hat{f}_nT_n^\\prime(x)\n", "$$\n", "\n", "**4)** Show using $\\theta$ variable that:\n", "$$\n", " \\mathbb{I}_N[f]^\\prime(\\cos(\\theta)) = \\frac1{\\sqrt{1-\\cos(\\theta)^2}}\\sum_{n=0}^Nn~\\hat{f}_n\\sin(n\\theta) \\text{ with } \\theta \\in ]0,\\pi[\n", "$$\n", "Then show that:\n", "$$\n", " \\mathbb{I}_N[f]^\\prime(1) = \\sum_{n=0}^Nn^2~\\hat{f}_n\\\\\n", " \\mathbb{I}_N[f]^\\prime(-1) = \\sum_{n=0}^N(-1)^{n+1}n^2~\\hat{f}_n\n", "$$\n", "\n", "**5)** Write a function which computes the values $\\{\\mathbb{I}_N[f]^\\prime(x_j)\\}_j$ given the coefficients of $\\mathbb{I}_N[f]^\\prime$. Complexity should be a $\\mathcal{O}(N\\log(N))$.\n", "\n", "To save some time, we can use the following function which relies on the [Discrete Sine Transform](https://en.wikipedia.org/wiki/Discrete_sine_transform) and its FFT algorithm ($\\mathcal{O}(N\\log(N))$ complexity) and computes $M$ values $\\{v(\\theta_j)\\}_j$ for a given set of $M$ coefficients $\\{\\hat{v}_m\\}_m$, related to each other through the following formula:\n", "$$\n", " v(\\theta_j) = \\sum_{m=1}^M\\hat{v}_m\\sin(m\\theta_j) \\text{ with } \\theta_j = \\frac{\\pi (j + 1)}{M + 1} \\text{ and } j = 0,...,M - 1\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy as scp\n", "import scipy.fftpack\n", "\n", "def InvDST (x):\n", " return scp.fftpack.idst (x / 2, type = 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**6)** Compare by plotting the exact value of $f^\\prime$ and the approximation $\\mathbb{I}_N[f]^\\prime$ at the grid points in the following cases:\n", "* $f(x) = \\cos((x + 1)\\pi) + \\sin(2(x + 1)\\pi)$\n", "* $f(x) = \\begin{cases}\n", " \\frac{\\pi}{4} & \\text{if } -1 \\leq x < 0 \\\\\n", " -\\frac{\\pi}{4} & \\text{if } 0 \\leq x < 1\n", " \\end{cases}$\n", "* $f(x) = 1/(1+25x^2)$\n", "\n", "If we have the time we can also plot the continuous function $\\mathbb{I}_N[f]^\\prime$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "\n", "\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }