{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TD6: Solving the biharmonic equation\n",
"\n",
"In this TD we want to solve a **biharmonic equation** in a 2D spatial domain $\\Omega = \\mathbb{R}\\times[-1, 1]$. We impose Dirichlet and Neumann boundary conditions and we will assume a $2\\pi$-periodicity of the solution along $x$-axis.\n",
"\n",
"For practical purposes, we will first be interested in the following $x$-invariant problem:\n",
"$$\n",
" \\begin{cases}\n",
" \\Delta^2 f (x, y) &= 0 \\text{ with } (x, y) \\in \\Omega\\\\\n",
" f (x, \\pm 1) &= 0 \\text{ with } x \\in [0,2\\pi[\\\\\n",
" \\partial_yf (x, \\pm 1) &= 1 \\text{ with } x \\in [0,2\\pi[\n",
" \\end{cases}\n",
"$$\n",
"whose exact solution is $f(x,y) = y\\left(y^2-1\\right)/2$.\n",
"\n",
"We will adopt a **spectral method** to solve this problem, approximating the solution by an expansion over Chebyshev polynomials for $y$-variable and over Fourier series for $x$-variable:\n",
"$$\n",
" f^{\\text{appr}}(x, y)= \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}\\hat{f}_{k,n}T_n(y)e^{ikx}\n",
"$$\n",
"The goal is to find $\\{\\hat{f}_{k,n}\\}$ spectral coefficients so that $f^{\\text{appr}}$ can be a solution of the biharmonic equation **and** so that $f^{\\text{appr}}$ can satisfy the boundary conditions. Unfortunately, in a general case, this goal is a pipe dream as the solution has no reason to be of the form of $f^{\\text{appr}}$, eventually we will need to make some compromises...\n",
"\n",
"Let us first define the orders of the expansion $N_x$ and $N_y$ with Python:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"Nx = 30\n",
"Ny = 30"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 1: A first natural approach: solving the 4^{th} order equation\n",
"\n",
"In 2D cartesian coordinates, the bilaplacian operator takes the following form:\n",
"$$\n",
" \\Delta^2 = \\frac{\\partial^4}{\\partial x^4} + 2\\frac{\\partial^4}{\\partial x^2\\partial y^2} + \\frac{\\partial^4}{\\partial y^4}\n",
"$$\n",
"Let us write the action of this operator on $f^{\\text{appr}}$:\n",
"$$\n",
" \\Delta^2f^{\\text{appr}} (x, y) = \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}\\left(k^4\\hat{f}_{k,n}T_n(y) - 2k^2\\hat{f}_{k, n}\\frac{\\partial^2}{\\partial y^2}T_n (y) + \\hat{f}_{k, n}\\frac{\\partial^4}{\\partial y^4}T_n (y)\\right)e^{ikx}\n",
"$$\n",
"In this case, the 2D biharmonic equation reduces to a set of $N_x$ independent 1D 4^{th} order equations. For a given Fourier mode $k$, we thus need to solve:\n",
"$$\n",
" \\sum_{n=0}^{N_y}k^4\\hat{f}_{k,n}T_n(y) - 2k^2\\hat{f}_{k, n}\\frac{\\partial^2}{\\partial y^2}T_n (y) + \\hat{f}_{k, n}\\frac{\\partial^4}{\\partial y^4}T_n (y) = 0\n",
"$$\n",
"which can be expressed as a matrix inversion problem:\n",
"$$\n",
" A_kX_k = B_k\n",
"$$\n",
"with $X_k = \\{\\hat{f}_{k, n}\\}_{n\\in\\{0, ..., N_y\\}}$, $A_k$ a matrix of shape $(N_y + 1)\\times(N_y + 1)$ which encodes the action of the previous 4^{th} order operator (in spectral space) as well as the boundary conditions through **tau method**, and $B_k$ a vector of length $N_y + 1$ whose non-zero coefficients encode the spectral values of the solution at the boundaries.\n",
"\n",
"For convenience, we recall that one can easily obtain the set of $\\{k\\}$ wave numbers using the following Python code:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kVals = 2 * np.pi * np.fft.fftfreq (Nx, 2 * np.pi / Nx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also expect the $\\{\\hat{f}_{k,n}\\}$ coefficients to be defined in Python as a matrix of complex numbers of shape $N_x\\times(N_y+1)$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fkn = np.zeros ((Nx, Ny + 1), dtype = complex)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Matrix problem conversion\n",
"\n",
"We first want to build the set of $\\{A_k\\}$ matrix. In the following, $k$ will refer to a given Fourier mode.\n",
"\n",
"**1)** Using *ChebDiffMatrix* function below which takes the order of Chebyshev expansion $N_y$ as a parameter and returns the associated Chebyshev differentiation matrix in spectral space (of shape $(N_y + 1)\\times(N_y+1)$), build a first matrix $\\tilde{A}_k$ which encodes the following differential operator in Chebyshev space:\n",
"$$\\Delta^2_k = k^4\\text{Id} - 2k^2\\frac{\\partial^2}{\\partial y^2} + \\frac{\\partial^4}{\\partial y^4}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def ChebDiff (x):\n",
" ret = np.zeros (len (x))\n",
" for k in range (len (x) - 1):\n",
" for p in range (k + 1, len (x)):\n",
" if (p + k) % 2 == 1:\n",
" ret [k] = ret [k] + 2 * p * x [p]\n",
" ret [0] = ret [0] / 2\n",
" return ret\n",
"\n",
"def ChebDiffMatrix (Ny):\n",
" D = np.zeros ((Ny + 1, Ny + 1))\n",
" for i in range (Ny + 1):\n",
" Xi = np.zeros (Ny + 1)\n",
" Xi [i] = 1.0\n",
" D [:, i] = ChebDiff (Xi)\n",
" return D\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**2)** To enforce boundary conditions, we will use **tau method** which fixes the values of the numerical solution at boundaries at the expense of higher order linear relations of the previous differential operator $\\Delta^2_k$. Using the following relations:\n",
"$$\\begin{cases}\n",
" f^{\\text{appr}}(x, 1) &= \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}\\hat{f}_{k,n}e^{ikx}\\\\\n",
" f^{\\text{appr}}(x, -1) &= \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}(-1)^n\\hat{f}_{k,n}e^{ikx}\\\\\n",
" \\partial_yf^{\\text{appr}}(x, 1) &= \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}n^2\\hat{f}_{k,n}e^{ikx}\\\\\n",
" \\partial_yf^{\\text{appr}}(x, -1)& = \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}(-1)^{n+1}n^2\\hat{f}_{k,n}e^{ikx}\n",
" \\end{cases}\n",
"$$\n",
"change $\\tilde{A}_k$ last four rows (which involve higher order linear relations of the previous differential operator $\\Delta^2_k$), and define the relevant $B_k$ vector so that $f^{\\text{appr}}$ satisfies the boundary conditions through the linear systems: $\\forall k, A_kX_k = B_k$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Matrix inversion and plotting\n",
"\n",
"Now we have properly formulated the initial problem into a matrix problem, we can solve it and plot the result.\n",
"\n",
"**3)** Invert the matrix problem using for example [scipy.linalg.lu_factor](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html) and [scipy.linalg.lu_solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html) routines to get $\\{\\hat{f}_{k,n}\\}$ coefficients."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import scipy.linalg as la\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**4)** Using *fEval* function which is defined below and computes for a given set of $\\{\\hat{f}_{k,n}\\}$ coefficients the values of $f^{\\text{appr}}$ at given points defined by their $x$ and $y$ coordinates, plot the solution using [pcolor](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.pcolor.html) routine. Then compare it to the exact solution $f(x,y) = y\\left(y^2-1\\right)/2$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def ChebyshevInterpolant (fnVals, x):\n",
" return cheb.chebval (x, fnVals)\n",
"\n",
"# This function expects fkn to be a matrix of coefficients whose first dimension refers to the Fourier index k\n",
"# and whose second dimension refers to the Chebyshev index n.\n",
"def fEval (fkn, x, y):\n",
" ret = np.zeros (x.shape, dtype = complex)\n",
" k = 2 * np.pi * np.fft.fftfreq (len (fkn [:, 0]), 2 * np.pi / len (fkn [:, 0]))\n",
" for i in range (fkn.shape [0]):\n",
" ret = ret + ChebyshevInterpolant (fkn [i, :], y) * np.exp (1j * k [i] * x)\n",
" return ret\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 2: Another approach: two nested 2^{nd} order problems\n",
"\n",
"The previous method can show some unstabilities on some more sophisticated problems. In this part we will investigate another approach which relies on the decomposition of the initial fourth order problem in **two second order nested problems**:\n",
"$$\n",
"\\begin{cases}\n",
" \\Delta f(x,y) &= g (x,y) &\\text{ with } (x, y) \\in \\Omega\\\\\n",
" \\Delta g (x,y) &= 0 &\\text{ with } (x, y) \\in \\Omega\\\\\n",
" f (x, \\pm 1) &= 0 &\\text{ with } x \\in [0,2\\pi[\\\\\n",
" \\partial_yf (x, \\pm 1) &= 1 &\\text{ with } x \\in [0,2\\pi[\n",
"\\end{cases}\n",
"$$\n",
"To get $f$, we need to know $g$ as $f$ is involved in a Poisson problem where $g$ plays the role of the source term. The issue is we **do not have any boundary conditions on $g$**. Thus, we have no choice but to find a general solution of $\\Delta g = 0$ for arbitrary boundary conditions, and then build $f$ from this general solution so $f$ can satisfy the boundary conditions.\n",
"\n",
"### Finding a general solution of $\\Delta g = 0$ on $\\Omega$\n",
"\n",
"$\\Omega$ boundaries are $\\partial\\Omega = \\mathbb{R}\\times\\{1\\}\\cup\\mathbb{R}\\times\\{-1\\}$, but from a numerical point of view, as we assume $2\\pi$-periodicity along $x$-axis, the boundaries are a discrete number of points depending on the order of the Fourier expansion $N_x$. **Numerical boundaries are** actually:\n",
"$$\n",
" \\partial\\Omega_{\\text{num}} = \\{(\\frac{2\\pi m}{N_x}, 1), m = 0, ..., N_x - 1\\}\\cup\\{(\\frac{2\\pi m}{N_x}, -1), m = 0, ..., N_x - 1\\}\n",
"$$\n",
"Therefore, there are $2N_x$ boundary points from a numerical point of view. In the following, we will note $\\partial\\Omega_{\\text{num}} = \\{\\mathbf{r_i}\\}_{i\\in\\{0,..., 2N_x - 1\\}}$.\n",
"\n",
"To solve $\\Delta g = 0$ for an arbitrary set of Dirichlet boundary conditions, we need to solve $2N_x$ Laplace problems of the following form:\n",
"$$\n",
"\\begin{cases}\n",
" \\Delta g_i = 0 \\text{ on } \\Omega\\\\\n",
" g_i(\\mathbf{r}) = \\delta\\left(\\mathbf{r}-\\mathbf{r_i}\\right)\n",
"\\end{cases}\n",
"$$\n",
"where $\\delta$ is the 2D Dirac function, and $i\\in\\{0, ..., N_x-1\\}$. $g_i$ functions are called **Green's functions**. The linearity of the Laplace equation will easily lead to a general solution as a linear combination of the Green's functions.\n",
"\n",
"**5)** Exactly on the same model as in first part, solve these $2N_x-1$ Laplacian problems looking for a numerical solution of the following form:\n",
"$$\n",
" g_i^{\\text{appr}}(x, y)= \\sum_{k = \\left \\lfloor{-N_x/2}\\right \\rfloor + 1}^{\\left \\lfloor{N_x/2}\\right \\rfloor}\\sum_{n=0}^{N_y}\\hat{g}_{k,n}^{(i)}T_n(y)e^{ikx}\n",
"$$\n",
"We will make use of the *FourierTransform* routine which is defined below and computes for a given set $\\{h(\\frac{2\\pi m}{N_x})\\}_{m\\in\\{0, ..., N_x -1\\}}$ of $h$ values, the Fourier coefficients $\\{\\hat{h}_k\\}_k$ of the associated Fourier interpolant."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def FourierTransform (fVals):\n",
" return np.fft.fft (fVals) / len (fVals)\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving $\\Delta f= g$ on $\\Omega$\n",
"\n",
"Now that we can build a general numerical solution of $\\Delta g = 0$ equation, we can use it to solve the initial problem.\n",
"\n",
"Using our previous $g_i$ solutions, we can now solve the $2N_x$ following Poisson problems:\n",
"$$\n",
"\\begin{cases}\n",
" \\Delta f_i = g_i\\text{ on } \\Omega\\\\\n",
" f_i(\\mathbf{r}) = 0 \\text{ on } \\partial\\Omega\n",
"\\end{cases}\n",
"$$\n",
"Once it is done, we will be able to compute the derivatives of $\\{f_i\\}$ functions, and find the suitable linear combination of $\\{f_i\\}$ functions so that we can build a numerical solution $f$ with $\\partial_yf (x, \\pm 1) = 1$.\n",
"\n",
"**6)** Solve these problems on the same model as in first part."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**7)** For each of these solutions, compute the associated numerical derivatives $\\partial_yf_i$. We will need *ChebDiff* routine (given in question **1)**) which computes the derivative of a Chebyshev interpolant (in spectral space). This routine takes a set of $\\{\\hat{h}_n\\}$ Chebyshev coefficients and returns the set $\\{\\hat{h}^\\prime_n\\}$ of Chebyshev coefficients of the derivative of the associated interpolant."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Building the solution $f$ which meets the boundary conditions\n",
"\n",
"Now we need to find the coefficients $\\{\\alpha_i\\}_{i\\in\\{0,...,2N_x-1\\}}$ so that:\n",
"$$\n",
" \\begin{cases}\n",
" f = \\sum_{i = 0}^{2N_x-1}\\alpha_if_i\\\\\n",
" \\forall i\\in\\{0, ..., 2N_x-1\\}, f(\\mathbf{r_i}) = 0 \\text{ and } \\partial_yf(\\mathbf{r_i})=1\n",
" \\end{cases}\n",
"$$\n",
"\n",
"**8)** Build a **capacitance matrix** $C$ (also called **influence matrix**) of shape $2N_x\\times2N_x$ with:\n",
"$$\n",
" C_{ij} = \\partial_yf_j(\\mathbf{r_i})\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**9)** Find the $\\{\\alpha_i\\}_{i\\in\\{0,...,2N_x-1\\}}$ coefficients that satisfy the previous relations."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**10)** Plot the numerical solution $f$ you have just built and compare it to the exact solution."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bonus question: a \"true\" 2D problem\n",
"\n",
"Our initial problem was not a really 2D problem as the boundary conditions were $x$-invariant. But what we have done before should solve biharmonic problems with more sophisticated boundary conditions.\n",
"\n",
"**BONUS: 11)** Solve using the method of part 1 the following problem:\n",
"$$\n",
"\\begin{cases}\n",
" \\Delta ^2 f = 0 \\text{ on } \\Omega\\\\\n",
" f(x,1) = e\\left(x\\cos(x)+\\sin(x)\\right)\\\\\n",
" f(x, -1) = e^{-1}\\left(x\\cos(x)-\\sin(x)\\right)\\\\\n",
" \\partial_yf(x,1) = e\\left(x\\cos(x)+2\\sin(x)\\right)\\\\\n",
" \\partial_yf(x,-1) = e^{-1}x\\cos(x)\n",
"\\end{cases}\n",
"$$\n",
"Compare it to the exact solution: $f(x,y) = x\\cos(x)e^y + y\\sin(x)e^y$. Is this solution $2\\pi$-periodic? Can you see a Gibbs phenomenon in your numerical solution?"
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}