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