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