Five_point_Laplacian.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.sparse\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The 5-point Laplacian\n",
    "The function below sets up a matrix that approximates the Laplacian operator in two dimensions:\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla^2 u = u_{xx} + u_{yy}.\n",
    "\\end{align}\n",
    "\n",
    "The formula for the finite difference approximation is\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla^2 u(x_i,y_j) \\approx \\frac{U_{i+1,j} + U_{i-1,j} + U_{i,j+1} + U_{i,j-1} - 4 U_{ij}}{h^2} = AU.\n",
    "\\end{align}\n",
    "\n",
    "Here we construct the matrix $A$, assuming that $U$ is order row-wise.  The constant $h$ is the grid spacing (which is assumed to be the same in the $x$ and $y$ directions).  The domain is assumed to be the square $[a,b]^2$ and it is discretized using $m$ points (not counting boundary values) in each direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def five_pt_laplacian_sparse(m,a,b):\n",
    "    \"\"\"Construct a sparse finite difference matrix that approximates the Laplacian.\"\"\"\n",
    "    e=np.ones(m**2)\n",
    "    e2=([1]*(m-1)+[0])*m\n",
    "    e3=([0]+[1]*(m-1))*m\n",
    "    h=(b-a)/(m+1)\n",
    "    A=scipy.sparse.spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)\n",
    "    A/=h**2\n",
    "    A = A.todia()\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, this matrix has the following form:\n",
    "\n",
    "\\begin{align}\n",
    "  A = \\frac{1}{h^2}\\begin{bmatrix} T & I & & & \\\\ I & T & I & & \\\\ & I & T & I & \\\\ & & \\ddots & \\ddots & \\ddots \\\\ & & & I & T\n",
    "  \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "where $I$ is the $m \\times m$ identity matrix and $T$ is a tridiagonal $m\\times m$ matrix:\n",
    "\n",
    "\\begin{align}\n",
    "  T = \\begin{bmatrix} -4 & 1 & & & \\\\ 1 & -4 & 1 & & \\\\ & 1 & -4 & 1 & \\\\ & & \\ddots & \\ddots & \\ddots \\\\ & & & 1 & -4\n",
    "  \\end{bmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This structure can be seen using the `spy` function, which shows where the nonzero entries are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x108e50050>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD7CAYAAABOrvnfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADrtJREFUeJzt3WusZXV5x/HvU8DWKSmENI5aJhwuxVhiiZQaQmuYWkim\nlIIvDEI0wpD4qrZoDBQwqedlK1GxbXihXApR6AUJoQmmDNoxNqQEcBimDBQRRxkMM0YEYxMFOk9f\n7DUnZ47nus76r33W+n8/yUn2Ze39/Neeec66nPXb/8hMJNXhV6Y9AEn9seGlitjwUkVseKkiNrxU\nERteqkjvDR8R2yLimYj4TkT8VeFaWyLiPyLiqYj474j4y5L1mppHRcSuiPi3HmodHxH3RMTTEbE3\nIs4pXO/65rPcExF3RcSvdvz+t0XEgYjYM++xEyJiR0Q8GxEPRsTxhevd2HyeuyPi3og4rmS9ec99\nMiIORcQJXdVbTK8NHxFHAf8AbAN+B7g8It5ZsOTrwCcy8wzgHODPC9cDuBrYC/RxgcMXgAcy853A\n7wJPlyoUETPAR4GzMvNdwFHAZR2XuZ3J/435rgN2ZObpwNeb+yXrPQickZlnAs8C1xeuR0RsAS4A\nvt9hrUX1vYV/D/BcZu7LzNeBfwIuKVUsM1/KzCea2z9j0hBvL1UvIk4ELgRuAaJUnabWccB7M/M2\ngMx8IzNfLVjyp0x+gW6KiKOBTcCLXRbIzG8BP1nw8MXAHc3tO4D3l6yXmTsy81Bz9xHgxJL1Gp8D\nru2qznL6bvjfAl6Yd39/81hxzRbq3Uz+EUv5PHANcGilBTtwMvCjiLg9Ir4dEV+KiE2limXmy8Bn\ngR8APwReycyHStWbZ3NmHmhuHwA291DzsKuAB0oWiIhLgP2Z+WTJOof13fBTuY43Io4F7gGubrb0\nJWpcBBzMzF0U3ro3jgbOAm7OzLOA/6Xb3d0jRMSpwMeBGSZ7ScdGxIdK1VtMTq4D7+X/UER8Cngt\nM+8qWGMTcAPw6fkPl6oH/Tf8i8CWefe3MNnKFxMRxwBfBb6cmfcVLHUucHFEfA+4G3hfRNxZsN5+\nJluGR5v79zD5BVDK2cDDmfnjzHwDuJfJOpd2ICLeChARbwMOli4YEVcyOTQr/QvtVCa/QHc3/29O\nBB6PiLeUKth3wz8G/HZEzETEm4APAveXKhYRAdwK7M3Mm0rVAcjMGzJzS2aezORk1jcy8yMF670E\nvBARpzcPnQ88Vaoe8AxwTkS8uflcz2dycrK0+4ErmttXACV/aRMR25gcll2SmT8vWSsz92Tm5sw8\nufl/s5/JSdFyv9Qys9cf4E+A/wGeA64vXOsPmRxPPwHsan629bCO5wH391DnTOBRYDeTLe5xhetd\ny+SXyh4mJ9CO6fj972ZyfuA1Jud6tgMnAA8xOWP+IHB8wXpXAd9hcrb88P+XmwvU+8Xh9Vvw/PPA\nCSX/DaMpJKkCXmknVcSGlypiw0sVseGlihxd4k0jwjOB0hRl5qIX8LiFl2pS6O+1SXMl5BLPz5b8\nW6P1rDeUel3XWqn3Wu/SN1ck3cQkJnlLZv7t4sud/Y9wysyRj542E3Hp1iMfe35f5mNXLv+6xRz5\nur61XT9pGlo1/Lxc+/lMro9/NCLuz8xF8tinzMC/nHfkY7PA7ElHPnbpKl63mIWv61vb9ZP61/YY\nfp259q0ty7a2s99yW/st1/v6WW+gtVo3/Dpz7Vtblm0nM3f2WnDk62e9YdaC9g3vn92kAWp70m5V\nufaImIXTZibHtFuZwq6uNHoRsZVVNlfbhp/LtTOJ+30QuHzhQpk5OzlbvfAElqSuNIcFOwEi4tPL\nLduq4TPzjYj4GPDvTP4sd+viZ+glbSSt/w6fmV8DvtbhWCQV5qW1UkWKhGeO9Py+1V108vy+bl7X\nt6GMU6LMV1wdTsvlEokdSWWs1Hvu0ksVseGlihQ/hm+beuv7dX0byjg1Lj2ctGubeuv7dX0byjg1\nJu7SSxWx4aWK2PBSRWx4qSI2vFQRG16qiA0vVcSGlyoywrTcUNJrQxmnxsS0nDQipuUkzelhl74d\nQzeLG8o4tTGtZ265LcCdwFuYfE/9FzPz77oamKGbpQxlnNqI1rOFfx34RGY+ERHHAo9HxA6/vVba\nuFofw2fmS5n5RHP7Z8DTwNu7Gpik7nVy0q6ZkOLdwCNdvJ+kMtZ90q7Znb8HuLrZ0s9/bnbe3Z39\nT+oojV8fU00dLnQM8FXgy5l538LnM3N2Pe8vaWVrmWqq9S59RARwK7A3M29q+z6S+rOeY/g/AD4M\n/FFE7Gp+tnU0LkkFrGduuf/EK/WkQdmwV9oZulmu/hDGqY3I8Iw0IoZnJM2x4aWKbOBj+HbGnpZr\na+zrp9UZXcOPPy3X1tjXT6vhLr1UERteqogNL1XEhpcqYsNLFbHhpYrY8FJFbHipIiO88Gbsabm2\nxr5+Wg3TctKImJaTNMeGlypS/Bh+7CmtsafzhjJOrc56v6b6KOAxYH9m/tniS409pTX2dN5QxqnV\nWO8u/dXAXiaTSUra4NbzvfQnAhcCtwCejZcGYD1b+M8D1wCHOhqLpMJaHcNHxEXAwczc1cxrtdRy\ns3DaDMwymfpqyUUltdTH3HLnAhdHxIXArwG/ERF3ZuZH5i+UmbMRl26F2ZNa1pG0guJzy2XmDZm5\nJTNPBi4DvrGw2SVtPF1deONZemkA1n3hTWZ+E/hmB2ORVFgPabmxp7TGns4byji1GqblpBExLSdp\njg0vVWTDfuPNUFJafaflxv65qKwN2/DDSWn1nZYb++eiktyllypiw0sVseGlitjwUkVseKkiNrxU\nERteqogNL1VkA194M5SUVt9pubF/LirJtJw0IqblJM2x4aWKbOBj+HaGktIyZbe4oYxzqFo3fEQc\nz2TWmTOYfInlVZn5X10NrL2hpLRM2S1uKOMcpvVs4b8APJCZH4iIo4Ff72hMkgppO/PMccB7M/MK\ngMx8A3i1y4FJ6l7bk3YnAz+KiNsj4tsR8aWI2NTlwCR1r+0u/dHAWcDHMvPRiLgJuA746/kLTeaW\nm7OzmRJHUof6mFtuP7A/Mx9t7t/DpOGPkJmzLd9f0ir1MbfcS8ALEXF689D5wFNt3ktSf9Zzlv4v\ngK9ExJuA7wLbuxmSpFJaN3xm7gZ+v8OxSCpsdFfaDSelZcpu6fpDGOcwmZaTRsS0nKQ5NrxUkREe\nw7czlJTW2NNybY19/bpiw88ZSkpr7Gm5tsa+ft1wl16qiA0vVcSGlypiw0sVseGlitjwUkVseKki\nNrxUES+8mTOUlNbY03JtjX39umFaThoR03KS5tjwUkU8hl+nsae0xp7OG8o4u7KeueWuBz4MHAL2\nANsz8xddDWw4xp7SGns6byjj7EarXfqImAE+CpyVme8CjgIu625Ykkpou4X/KfA6sCki/g/YBLzY\n2agkFdF2IoqXgc8CPwB+CLySmQ91OTBJ3Ws7e+ypwMeBGSazxv5rRHwoM7+yYLnZeXedW04qoI+5\n5c4GHs7MHzcF7wXOBY5oeOeWk8orPrcc8AxwTkS8OSKCydxye1u+l6SetD2G3w3cCTwGPNk8/MWu\nBiWpjPXMLfcZ4DMdjkVSYV5pt25jT2mNPZ03lHF2w7ScNCKm5STNseGlingMPyVDSWn1nZYb++cy\nbTb81AwlpdV3Wm7sn8t0uUsvVcSGlypiw0sVseGlitjwUkVseKkiNrxUERteqogX3kzNUFJafafl\nxv65TJdpOWlETMtJmmPDSxXxGH5ghpLSMmX3yzbCGJdt+Ii4DfhT4GAzpRQRcQLwz8BJwD7g0sx8\npcTgtJihpLRM2f2y6Y9xpV3624FtCx67DtiRmacDX2/uSxqAZRs+M78F/GTBwxcDdzS37wDeX2Bc\nkgpoc9Juc2YeaG4fADZ3OB5JBa3rpF1m5uG/+y3GueWk8krPLXcgIt6amS9FxNuAg0st6NxyUnml\n55a7H7iiuX0FcF+L95A0Bcs2fETcDTwMvCMiXoiI7cDfABdExLPA+5r7kgZg2V36zLx8iafOLzAW\nSYV5pd3gDCWlZcpu8drTHaNpOWlETMtJmmPDSxXxGL4SGyGptRpjTsu1tZZ1W2kJG74a009qrc6Y\n03JtrWXdHl92CXfppYrY8FJFbHipIja8VBEbXqqIDS9VxIaXKmLDSxXxwptqTD+ptTpjTsu1taZ1\nW/YCHdNy0oiYlpM0x4aXKuIxvJY15hQa1JfOW7Hhl5hf7kbgIuA14LvA9sx8teRANS1jTqFBbem8\n1ezSLza/3IPAGZl5JvAscH3XA5PUvRUbfrH55TJzR2Yeau4+ApxYYGySOtbFSburgAc6eB9Jha3r\npF1EfAp4LTPvWuL52Xl3nVtOKqD03HKHi1wJXAj88VLLOLecVN5a5pZr1fARsQ24BjgvM3/e5j0k\n9W/FY/hF5pe7Cvh74FhgR0TsioibC49TUgdW3MIvMb/cbQXGIqkwr7TTCsacQoPa0nmm5aQRMS0n\naY4NL1XEY3gVMZQ0WZ9puY3wmdjwKmQoabI+03LT/0zcpZcqYsNLFbHhpYrY8FJFbHipIja8VBEb\nXqqIDS9VxAtvVMhQ0mR9puWm/5mYlpNGxLScpDk2vFSRZY/hF5tmat5znwRuBH4zM18uN0TVZCMk\nylajz7TcWl630hIrnbS7nckXVt45/8GI2AJcAHx/5UFIazH9RNnq9JmWW8vrHl92iWV36RebZqrx\nOeDalQcgaSNZ8zF8RFwC7M/MJwuMR1JBa/o7fERsAm5gsjs/93CnI5JUzFovvDkVmAF2RwRMZo19\nPCLek5kHFy7s3HJSH3Y2P7tmVlpyTQ2fmXuAzYfvR8T3gN9b6iy9c8tJfdja/OzdB8+dtNySyx7D\nz5tm6vRmmqntCxbp/jI9ScUsu4VfYpqp+c+f0u1wJJXklXZSRUzLaYOZfqJsdfpMy63pdcteoGNa\nThoR03KS5tjwUkU8htcoDCFlN4S0nDQQQ0jZbfC0nKRxseGlitjwUkVseKkiNrxUERteqshUGj4i\ntlrPetaDyRdX9GdaW/it1rOe9aDvhvfCG43EEFJ2i41x18zkm2oWLrfS65Z6/ymm5SRNh2k5SWW2\n8JI2JrfwUkVseKkiNrxUERteqogNL1Xk/wETjwUC0E67IwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108d61dd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = five_pt_laplacian_sparse(4,-1.,1.)\n",
    "plt.spy(A)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}