{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# TD2: Spatial Discretization with Spectral Methods I\n", "\n", "When facing a spatial differential equation that we want to solve numerically, with spatial discretization methods we want to compute all values $f(x_0), ..., f(x_N)$ of a discretized function $f$ simultaneously. Spatial discretization methods are therefore opposed to temporal discretization methods 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 yet in solving a spatial differential equation 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": [ "# To draw matplotlib plots within this notebook.\n", "%matplotlib inline\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: Fourier 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", "In this TD we only 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", "Note that this function is a continuous function, therefore it can be evaluated at any point in $[0,2\\pi]$. 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(y) / N* where *y* is the vector of values of the set $\\{f(x_m)\\}_m$ and $N$ the length of *y*. 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)** Is $\\mathbb{I}_N[f]$ always a real function?\n", "\n", "We will use and assume an **odd value of $N$** in the following." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**4)** 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) = \\sin (2x)\\exp\\left(\\cos\\left(8x\\right)\\right)$\n", " \n", "Can you observe the Gibbs phenomenon for one of these functions?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: Spectral differentiation with Fourier representation\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 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) = \\sin (2x)\\exp\\left(\\cos\\left(8x\\right)\\right)$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have been dealing so far with Fourier interpolation and Fourier differentiation only for odd values of N. Indeed, it is easier to define a real interpolant when symmetric modes are involved.\n", "\n", "Dealing with even values of $N$ can be more technical, curious readers could take a look at [Steven G. Johnson notes](https://math.mit.edu/~stevenj/fft-deriv.pdf) to get a deeper understanding of what is behind the definition of a Fourier interpolant." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 3: Aliasing around an example...\n", "\n", "Aliasing effects are caused by the fact that on a discrete grid $\\{x_m = \\frac{2\\pi m}{N}, m = 0, ..., N - 1\\}$, all modes $\\{\\phi_{k + jN}(x) = e^{i(k+jN)x}\\}_{j\\in\\mathbb{Z}}$ are indistinguishable (they are *aliases* of one another). Let us see on an example what it can look like...\n", "\n", "Let us consider the following functions of $L^2([0,2\\pi])$:\n", "* $u(x) = \\cos (2x)$\n", "* $v(x) = \\sin (x)$\n", "* $w (x) = u (x)\\times v (x) = \\cos (2x) \\sin (x)$\n", "\n", "**1)** What is the Fourier series associated with $w$? What is the minimum odd value of $N$ so that $\\mathbb{I}_N[w] = w$?\n", "\n", "**2)** Compute $\\mathbb{I}_N[w]$ and then plot $u, v, w$ and their interpolants $\\mathbb{I}_N[u], \\mathbb{I}_N[v], \\mathbb{I}_N[w]$ on the same figures for $N = 7$, $N = 5$ and $N = 3$. Discuss what you observe." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**3)** Generally with spectral methods it is easier to compute products in physical space to avoid some convolutions in spectral space. In this case we talk about **pseudo-spectral methods**. On a numerical problem, what should we do generally before and after computing products $u \\times v$ in physical space?\n", "\n", "**BONUS: 4)** Propose a *de-aliasing* algorithm that computes $\\hat{w}_k$ coefficients given the sets of $\\hat{u}_k$ and $\\hat{v}_k$ coefficients and which avoids unwanted aliasing effects." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "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.7.7" } }, "nbformat": 4, "nbformat_minor": 2 }