A partial differential equation (PDE) of the function $u(x_1,\\dots,x_n)$ is an equation of the form:\n", "\\begin{equation}\n", " F\\big( \\underbrace{x_1,\\dots,x_n}_{\\text{variables}} , \\underbrace{u}_{\\text{function}}, \\underbrace{ \\frac{\\partial u}{\\partial x_1}, \\dots, \\frac{\\partial u}{\\partial x_n}}_{\\text{first derivatives}}, \\underbrace{ \\frac{\\partial^2 u}{\\partial x_1^2}, \\frac{\\partial^2 u}{\\partial x_1\\partial x_2}, \\dots }_{\\text{second derivatives}}, \\underbrace{\\dots}_{\\text{higher order derivatives}} \\big) = 0 \n", "\\end{equation}\n", "\n", "
We will focus on linear $2^{nd}$ order PDEs in 2D, which can be written in the form: \n", "\\begin{equation}\n", "a \\frac{\\partial^2 u}{\\partial x^2} + b \\frac{\\partial^2 u}{\\partial x \\partial y} + c \\frac{\\partial^2 u}{\\partial y^2} + \\underbrace{\\dots}_{\\text{lower terms}} = 0\n", "\\end{equation}\n", "Like for conic equations, we can define three classes of linear $2^{nd}$ order PDEs, based on the value of discriminant $b^2 - 4 a c $ : \n", "
The Poisson equation is a $2^{nd}$ order PDE of the elliptic type: \n", "\\begin{equation}\n", "\\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2} = f(x,y)\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a 2D electrostatic problem the Poisson equation can be written as\n", "\\begin{equation}\n", "\\nabla^2 V(x,y) = - \\rho( x, y )\n", "\\end{equation}\n", "where $\\nabla^2=\\sum_i\\frac{\\partial^2}{\\partial x_i^2}$ is the Laplace operator, $V$ is the potential, and $\\rho$ is the charge density. We note that the Poisson equation differs from the Lapace equation because of the presence of the source term $\\rho$, which makes the equation inhomogeneous. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
For the aforementioned source term, it is easy to verify that the solution of the 2D Poisson equation is:\n", "\\begin{equation}\n", "V(x,y)=\\frac{L^2}{40\\pi^2} sin\\left( 6 \\pi \\frac{x}{L} \\right) sin\\left( 2 \\pi \\frac{y}{L} \\right)\n", "\\end{equation}\n", "Let's plot $\\rho(x,y)$ and $V(x,y)$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define the periodicity\n", "L = 10.0\n", "# Define the number of points used to discretize the source and the potential on a NxN mesh\n", "N = 25\n", "\n", "# Create the list that contains the x positions \n", "x = []\n", "for i in range(N) : \n", " x.append(float(i)/float(N)*L)\n", "\n", "# Create the list that contains the y positions \n", "y = []\n", "for j in range(N) : \n", " y.append(float(j)/float(N)*L)\n", "\n", "# Show the arrays\n", "print(\"x = \", x)\n", "print(\"y = \", y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Access elements of the array x (Python experts can skip)\n", "print(\"x[0] = \", x[0]) \n", "print(\"x[1] = \", x[1]) \n", "print(\"x[N-1] = \", x[N-1]) \n", "print(\"x[-1] = \", x[-1]) \n", "print(\"x[-N] = \", x[-N]) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# create X and Y grids to easily handle gridpoints\n", "import numpy as np\n", "X, Y = np.meshgrid(x, y)\n", "print(\"X = \", X)\n", "print(\"Y = \", Y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define rho and V on the grid\n", "import numpy as np\n", "\n", "# rho \n", "rho = np.empty([N,N]) \n", "rho = np.sin(X/L*np.pi*6.0)*np.sin(Y/L*np.pi*2.0)\n", "\n", "# V analytical solution \n", "Van = np.empty([N,N]) \n", "Van = (np.sin(X/L*np.pi*6.0)*np.sin(Y/L*np.pi*2.0))/(40.0*np.pi*np.pi/L/L)\n", "\n", "print(\"rho = \", rho)\n", "print(\"Van = \", Van)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# method to plot a 3D function \n", "def plot3D( X, Y, Z ) :\n", " from mpl_toolkits.mplot3d import Axes3D\n", " import matplotlib.pyplot as plt\n", " from matplotlib import cm\n", " from matplotlib.ticker import LinearLocator, FormatStrFormatter\n", " import numpy as np\n", "\n", "\n", " fig = plt.figure()\n", " ax = fig.gca(projection='3d')\n", " \n", " # Plot the surface.\n", " surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,\n", " linewidth=0, antialiased=False)\n", "\n", " # Customize the z axis.\n", " ax.set_zlim(-1.01, 1.01)\n", " ax.zaxis.set_major_locator(LinearLocator(10))\n", " ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))\n", "\n", " # Set labels\n", " ax.set_xlabel('x')\n", " ax.set_ylabel('y')\n", " ax.set_zlabel('z')\n", "\n", " # Add a color bar which maps values to colors.\n", " fig.colorbar(surf, shrink=0.5, aspect=5)\n", "\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot3D( X, Y, rho)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot3D( X, Y, Van)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Both the potential and the source are 2D periodic functions, with periodicity $[L,L]$, and can be therefore expanded in plane waves using the following Fourier expansion:\n", "\\begin{equation}\n", " \\rho(\\mathbf{r}) = \\sum_\\mathbf{G} \\rho(\\mathbf{G}) e^{i\\mathbf{G}\\cdot\\mathbf{r}}\n", "\\end{equation}\n", "where the Fourier components $\\rho(\\mathbf{G})$ are obtained as:\n", "\\begin{equation}\n", " \\rho(\\mathbf{G}) = \\frac{1}{L^2}\\int d\\mathbf{r} \\rho(\\mathbf{r}) e^{-i\\mathbf{G}\\cdot\\mathbf{r}}\n", "\\end{equation}\n", "$\\mathbf{r}=(x,y)$, $\\mathbf{G}=\\frac{2\\pi}{L}(n,m)$, with $n=0,\\pm 1, \\pm 2, \\dots $ and $m=0,\\pm 1, \\pm 2, \\dots $." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Poisson equation in Fourier space becomes a simple multiplication\n", "\\begin{equation}\n", "G^2 V(\\mathbf{G}) = \\rho(\\mathbf{G})\n", "\\end{equation}\n", "where $G^2 = \\mathbf{G}\\cdot\\mathbf{G}=G_x^2+G_y^2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can therefore follow this protocol to obtain $V(x,y)$ : \n", "