{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TD3: Spatial Discretization with Spectral Methods II\n",
"\n",
"In this TD we continue to explore spectral methods to represent a given function and its derivatives. We will focus here on **Chebyshev representation** which most of the time can be seen as a reasonable alternative to Fourier representation in the case where the function that we want to approximate on a compact domain is regular but non-periodic.\n",
"\n",
"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: 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",
"### Chebyshev basis\n",
"\n",
"Chebyshev interpolation relies on Chebyshev polynomials. Let us first see what the Chebyshev polynomials look like.\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": [
"**1)** 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": [
"**2)** 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": [
"**3)** 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 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 cheb_transform(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": [
"**4)** 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": [
"**5)** 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 interpolants seem to be relevant approximations? Can you see a Gibbs phenomenon for one of these interpolants?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 2: Spectral differentiation with Chebyshev 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 Chebyshev representation of a given function $f$.\n",
"\n",
"### Chebyshev 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_{n=0}^N\\hat{f}_nT_n^\\prime(x)\n",
"$$\n",
"\n",
"Now we can wonder if this function is still a good approximation of $f^\\prime$.\n",
"\n",
"**1)** 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",
"**2)** Write a function which computes the values $\\{\\mathbb{I}_N[f]^\\prime(x_j)\\}_j$ given the coefficients of $\\mathbb{I}_N[f]$. 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": [
"def inv_dst(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": [
"**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 + 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",
"Discuss what you observe."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## [BONUS] Part 3: Lagrange interpolation and the Runge phenomenon\n",
"\n",
"Chebyshev interpolation can be seen as a polynomial interpolation of a function based on grid points that are not evenly spaced. Indeed these grid points are cosine projections of evenly spaced points. One can wonder why we could not use an evenly spaced grid to define a polynomial interpolation. Let us see what can happen.\n",
"\n",
"We want to approximate the function $f(x) = 1 / (1 + 25x^2)$ on $[-1,1]$ using a polynomial interpolation on a regular grid.\n",
"\n",
"For a given $N$ we define a grid of $N + 1$ points $\\{x_j\\}_{0\\leq j \\leq N}$ on $[-1, 1]$ by:\n",
"$$\n",
" x_j = -1 + \\frac{2j}{N} \\text{ with } j=0,...,N\n",
"$$\n",
"\n",
"A polynomial interpolation on this grid can be defined using Lagrange polynomials. We recall that the Lagrange polynomials $\\{\\mathcal{l}_k\\}_k$ associated to a given grid of $N + 1$ points $\\{x_j\\}_{0\\leq j \\leq N}$ are defined as:\n",
"$$\n",
" \\mathcal{l}_k (x) = \\prod_{0 \\leq j \\leq N, j\\neq k}\\frac{x - x_j}{x_k-x_j}\n",
"$$\n",
"Thus, they verify the relation: $\\mathcal{l}_k (x_j) = \\delta_{j,k}$\n",
"\n",
"The polynomial interpolant of a function $f$ on the grid $\\{x_j\\}_{0\\leq j \\leq N}$ is then simply:\n",
"$$\n",
" \\mathbb{I}_N[f](x) = \\sum_{k=0}^N\\hat{f}_k\\mathcal{l}_k(x) \\text{ with } \\hat{f}_k = f(x_k)\n",
"$$\n",
"\n",
"**BONUS: 1)** $\\mathbb{I}_N[f]$ can be written $\\mathbb{I}_N[f](X) = \\sum_{k=0}^Na_kX^k$. Using [scipy.interpolate.lagrange](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html) define a function that computes the $\\{a_k\\}_{0\\leq k \\leq N}$ coefficients associated to the function $f(x) = 1 / (1 + 25x^2)$ for a given value of $N$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**BONUS: 2)** Plot on $[-1,1]$ for several values of $N$ the function $f$ and its interpolant $\\mathbb{I}_N[f]$. What happens when you increase the order of interpolation $N$? Conclude on the interest of evenly spaced grids for polynomial interpolation."
]
},
{
"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
}