{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD4: 2D Laplace equation\n",
    "\n",
    "In this TD we want to solve the following [Dirichlet problem](https://en.wikipedia.org/wiki/Dirichlet_problem) on a bounded 2D spatial domain $\\Omega$:\n",
    "$$\n",
    "    \\begin{cases}\n",
    "        \\Delta f = 0 \\text{ on } \\Omega\\\\\n",
    "        f = g \\text{ on } \\partial\\Omega\n",
    "    \\end{cases}\n",
    "$$\n",
    "where $\\partial\\Omega$ is the frontier of the domain $\\Omega$ and $g$ is a given function which determines the solution at the boundaries of the domain. From a physical point of view, the Laplace equation $\\Delta f = 0$ is the steady-state heat equation, thus $f$ can be interpreted as a field of temperature and $g$ as the temperature imposed at the boundaries of the domain $\\Omega$.\n",
    "\n",
    "In the following we will work on a squared domain $\\Omega = ]-\\frac{L}{2}, \\frac{L}{2}[^2$ and we will use a **finite-difference method** to solve this equation. We define a five-point stencil noted $\\Delta_5$ and a nine-point stencil noted $\\Delta_9$ as follows:\n",
    "$$\n",
    "\\Delta_5 f (x, y) = \\frac1{h^2}\\left(f \\left(x + h, y\\right) + f \\left(x - h, y\\right) + f \\left(x, y +h\\right) + f\\left(x, y -h\\right) - 4f\\left(x, y\\right)\\right)\n",
    "$$\n",
    "$$\n",
    "\\Delta_9 f (x, y) = \\frac1{6h^2}\\left(4\\left(f \\left(x + h, y\\right) + f \\left(x - h, y\\right) + f \\left(x, y +h\\right) + f\\left(x, y -h\\right)\\right) + f \\left(x+h, y+h\\right) + f \\left(x-h, y+h\\right) + f \\left(x-h, y-h\\right) + f \\left(x+h, y-h\\right) - 20f\\left(x, y\\right)\\right)\n",
    "$$\n",
    "We can use these operators for a discrete function $f$ defined on a regular grid which discretizes the spatial domain $\\Omega$. We just need to make sure that the $x$-discretization step is the same as the $y$-discretization step. These two operators can also be schematically represented by the two following arrays of coefficients:\n",
    "$$\n",
    "    \\begin{array}{|c|c|c|}\n",
    "    \\hline\n",
    "    &1&\\\\\n",
    "    \\hline\n",
    "    1&-4&1\\\\\n",
    "    \\hline\n",
    "    &1&\\\\\n",
    "    \\hline\n",
    "    \\end{array}\n",
    "    \\begin{array}{|c|c|c|}\n",
    "    \\hline\n",
    "    1/6&2/3&1/6\\\\\n",
    "    \\hline\n",
    "    2/3&-10/3&2/3\\\\\n",
    "    \\hline\n",
    "    1/6&2/3&1/6\\\\\n",
    "    \\hline\n",
    "    \\end{array}\n",
    "$$\n",
    "Positions of the coefficients in the array refer to the neighbors of a given point in the grid. This given point is represented by the middle coefficient of the array.\n",
    "\n",
    "To solve this kind of linear spatial problem we need to proceed in three steps:\n",
    "1. Define a grid which covers $\\Omega\\cup\\partial\\Omega$ and is suitable for the chosen solving method.\n",
    "2. Build the matrix which constrains the solution on the domain.\n",
    "3. Find the inverse of this matrix to solve the problem.\n",
    "\n",
    "We will need the following Python packages for this problem:"
   ]
  },
  {
   "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\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**1)** Defining $h$ as the step of discretization along $x$ and $y$ axis, build a squared grid $\\mathcal{G}$ which covers $\\Omega\\cup\\partial\\Omega$. We can use the [numpy meshgrid](https://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html) function to get useful arrays of $x$ and $y$ coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining $N$ as the number of points in the grid along $x$ or $y$ axis, now we need to impose $N\\times N$ constraints on the numerical solution $f$. The linearity of the problem naturally leads to the definition of a matrix of constraints $A$.\n",
    "\n",
    "With $(x_0, y_0)$ the bottom left corner of the grid $\\mathcal{G}$, we should have:\n",
    "$$\n",
    "    \\mathcal{G} = \\{(x_i = x_0 + ih, y_j = y_0 + jh), i = 0, ..., N - 1 \\text{ and } j = 0, ..., N-1\\} \n",
    "$$\n",
    "\n",
    "To solve our Dirichlet problem, let us approximate the laplacian operator $\\Delta$ by the **five-point stencil** $\\Delta_5$ previously defined. Now we want to find a solution $f$ defined on the grid $\\mathcal{G}$ so that:\n",
    "1. $\\forall (x_i, y_j) \\in \\mathcal{G} \\cap \\Omega, \\Delta_5 f(x_i, y_j) = 0$\n",
    "2. $\\forall (x_i, y_j) \\in \\mathcal{G} \\backslash \\Omega, f(x_i, y_j) = g(x_i, y_j)$\n",
    "where $g$ is a given function.\n",
    "\n",
    "We will store the values of the solution $f$ over the grid in a vector $X$ of length $N \\times N$ where $X_{iN + j}$ coordinate will be the value of the solution $f$ at the grid point $(x_i, y_j)$.\n",
    "\n",
    "**2)** Build a matrix $A$ of shape $(N \\times N, N\\times N)$ so that the previous constraints on the solution $f$ can be expressed as a linear problem:\n",
    "$$\n",
    "    AX = B\n",
    "$$\n",
    "with $B$ a vector of length $N \\times N$ that only depends on $\\{g(x_i, y_j)\\}_{i,j}$ values. We can use in the following: $g(x,y) = \\cos (2\\pi(x + y)/L)$ to test our codes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** Using [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) functions, solve the previous linear system through a LU decomposition of $A$. **Beware of the $\\mathcal{O}(M^3)$ complexity of LU decomposition algorithm** where $M$ is the order of the input matrix (in our case, $M = N\\times N$), computation time can be long!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4)** Plot the solution of the Dirichlet problem on the domain $\\Omega$. We may want to use [matplotlib.pyplot.pcolor](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.pcolor.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5)** Solve the same problem approximating the laplacian operator $\\Delta$ with the **nine-point stencil** $\\Delta_9$."
   ]
  },
  {
   "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
}
