{
 "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
}