lessons/05_relax/01_2D_Laplace_Equation.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:6830d350704d9d3e2c977a40f8c35cfc2798de02af58c100c4846236b398dc13"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth.  Based on [CFDPython](https://github.com/barbagroup/CFDPython), (c)2013 L.A. Barba, also under CC-BY license."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Relax and hold steady"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You've reached the fifth module of the open course [**\"Practical Numerical Methods with Python\"**](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about), titled *Relax and hold steady: elliptic problems*. We hope that you are enjoying the ride of [#numericalmooc](https://twitter.com/hashtag/numericalmooc) so far!\n",
      "\n",
      "We introduced finite-difference methods for partial differential equations (PDEs) in the [second module](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/tree/master/lessons/02_spacetime/), and looked at convection problems in more depth in [module 3](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/tree/master/lessons/03_wave/). \n",
      "\n",
      "In [module 4](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/tree/master/lessons/04_spreadout/) we explored diffusion-dominated problems.  It's been quite a trip so far and now it's time to stop worrying about **time** and relax.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Laplace's Equation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here is Laplace's equation in 2D:\n",
      "\n",
      "\\begin{equation}\n",
      "\\frac{\\partial ^2 p}{\\partial x^2} + \\frac{\\partial ^2 p}{\\partial y^2} = 0\n",
      "\\end{equation}\n",
      "\n",
      "See anything familiar?  What should jump out at you right away is the *strong* similarity to the diffusion equation that was the focus of [module 4](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/tree/master/lessons/04_spreadout/), but what's different?\n",
      "\n",
      "The Laplace Equation does not have a time dependence \u2014 there is no $p^{n+1}$.  Instead of tracking a wave through time as in previous modules, the Laplace equation calculates the equilibrium state of a system under the supplied boundary conditions.  \n",
      "\n",
      "If you have taken coursework in Heat Transfer, you will recognize the Laplace Equation as the steady-state heat equation."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Laplace's equation has the features typical of diffusion phenomena. For this reason, it has to be discretized with *central differences*, so that the discretization is consistent with the physics we want to simulate. \n",
      "\n",
      "You should be able to write out a second-order central difference scheme by now.  \n",
      "\n",
      "\\begin{equation}\n",
      "\\frac{p_{i+1, j}^n - 2p_{i,j}^n + p_{i-1,j}^n}{\\Delta x^2} + \\frac{p_{i,j+1}^n - 2p_{i,j}^n + p_{i, j-1}^n}{\\Delta y^2} = 0\n",
      "\\end{equation}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can rearrange this equation and solve for $p_{i,j}^{n+1}$ to get\n",
      "\n",
      "\\begin{equation}\n",
      "p_{i,j}^{n+1} = \\frac{\\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\\Delta x^2 + \\Delta y^2)}\n",
      "\\end{equation}\n",
      "\n",
      "Since there is no time dependence in the Laplace Equation, instead of calculating where the system will be at some time $t$, we will iteratively solve for $p_{i,j}^{n}$ until it meets a condition that we specify.  The system will reach equilibrium only as the number of iterations tends to $\\infty$, but we can approximate the equilibrium state by iterating until the change between one iteration and the next is *very* small.  \n",
      "\n",
      "This method of iterative solution is known as the **Jacobi method**.\n",
      "\n",
      "Using second-order central-difference schemes in both directions is the most widely applied method for the Laplace operator. It is also known as the **five-point difference operator**, alluding to its stencil.\n",
      "\n",
      "<br>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<img src=\"./figures/laplace.svg\">"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Figure 1: Laplace five-point stencil"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Very often we deal with meshes where $\\Delta x = \\Delta y$ and we can use this to further simplify Equation $(3)$.  \n",
      "\n",
      "If $\\Delta x = \\Delta y$, Equation $(3)$ becomes\n",
      "\n",
      "\\begin{equation}\n",
      "p^{n+1}_{i,j} = \\frac{1}{4} \\left(p^{n}_{i,j-1} + p^n_{i,j+1} + p^{n}_{i-1,j} + p^n_{i+1,j} \\right)\n",
      "\\end{equation}\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Boundary conditions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to solve Laplace's equation numerically by assuming an initial state of $p=0$ everywhere and then add the following boundary conditions:\n",
      "\n",
      "\\begin{equation}\n",
      "  \\begin{gathered}\n",
      "p=0 \\text{ at } x=0\\\\\n",
      "\\frac{\\partial p}{\\partial x} = 0 \\text{ at } x = L\\\\\n",
      "p = 0 \\text{ at }y = 0 \\\\\n",
      "p = \\sin \\left(  \\frac{\\frac{3}{2}\\pi x}{L} \\right) \\text{ at } y = H\n",
      "  \\end{gathered}\n",
      "\\end{equation}\n",
      "\n",
      "\n",
      "Where $L=1$ and $H=1$ are the sizes of the domain in the $x$ and $y$ directions, respectively."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One of the defining features of elliptic PDEs is that they are \"driven\" by the boundary conditions.  The boundaries are set to specific values and the solution relaxes to join the boundaries together.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First let's import our usual smattering of libraries (plus a few new ones!)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import pyplot\n",
      "import numpy\n",
      "%matplotlib inline\n",
      "from matplotlib import rcParams\n",
      "rcParams['font.family'] = 'serif'\n",
      "rcParams['font.size'] = 16"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "When we want to visualize 2D data, we can often use `pyplot.imshow` as we did in the Reaction-Diffusion assignment, but often a 3D projection plot can provide a more intuitive view of a system.  \n",
      "\n",
      "For projection plots, we have to import the `Axes3D` library from matplotlib.  We'll also grab the `cm` package, which provides different colormaps for visualizing plots.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mpl_toolkits.mplot3d import Axes3D\n",
      "from matplotlib import cm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rather than repeatedly setting up a 3D plotting environment to visualize each result, we can define a function to save time later.  \n",
      "\n",
      "Many examples of 3D projection plots are available in the mplot3d section of the [matplotlib gallery](http://matplotlib.org/gallery.html#mplot3d).  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plot2D(x, y, p):\n",
      "    '''Creates 3D projection plot with appropriate limits and viewing angle\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    x: array of float\n",
      "        nodal coordinates in x\n",
      "    y: array of float\n",
      "        nodal coordinates in y\n",
      "    p: 2D array of float\n",
      "        calculated potential field\n",
      "    \n",
      "    '''\n",
      "    fig = pyplot.figure(figsize=(11,7), dpi=100)\n",
      "    ax = fig.gca(projection='3d')\n",
      "    X,Y = numpy.meshgrid(x,y)\n",
      "    surf = ax.plot_surface(X,Y,p[:], rstride=1, cstride=1, cmap=cm.coolwarm,\n",
      "            linewidth=0, antialiased=False)\n",
      "\n",
      "    ax.set_xlim(0,1)\n",
      "    ax.set_ylim(0,1)\n",
      "    ax.set_xlabel('$x$')\n",
      "    ax.set_ylabel('$y$')\n",
      "    ax.set_zlabel('$z$')\n",
      "    ax.view_init(30,45)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Analytical Solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Laplace equation with the boundary conditions listed in $(4)$ have an analytical solution, given by\n",
      "\n",
      "\\begin{equation}\n",
      "p(x,y) = \\frac{\\sinh \\left( \\frac{\\frac{3}{2} \\pi y}{L}\\right)}{\\sinh \\left(  \\frac{\\frac{3}{2} \\pi H}{L}\\right)} \\sin \\left( \\frac{\\frac{3}{2} \\pi x}{L} \\right)\n",
      "\\end{equation}\n",
      "\n",
      "where $L$ and $H$ are the length of the domain in the $x$ and $y$ directions, respectively."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We used `numpy.meshgrid` to plot our 2D solutions to the heat equation in Module 4.  Here we again use it as a plotting aid.  `linspace` creates 1-row arrays which work well for defining $x$ and $y$ axes in line plots, but now we want the analytical solution at every point in our domain.  To do this, we create a meshgrid and then use the X and Y values from `numpy.meshgrid` in the analytical solution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def p_analytical(x, y):\n",
      "    X, Y = numpy.meshgrid(x,y)\n",
      "    \n",
      "    p_an = numpy.sinh(1.5*numpy.pi*Y / x[-1]) /\\\n",
      "    (numpy.sinh(1.5*numpy.pi*y[-1]/x[-1]))*numpy.sin(1.5*numpy.pi*X/x[-1])\n",
      "    \n",
      "    return p_an"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ok, let's try out the analytical solution and use it to test the `plot2D` function we wrote above.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx = 41\n",
      "ny = 41\n",
      "\n",
      "x = numpy.linspace(0,1,nx)\n",
      "y = numpy.linspace(0,1,ny)\n",
      "\n",
      "p_an = p_analytical(x,y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot2D(x,y,p_an)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It worked!  This is what the solution *should* look like when we're 'done' relaxing.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "How long do we iterate?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We noted above that there is no time dependence in the Laplace equation.  Given that, using a `for` loop with `nt` iterations doesn't make a lot of sense.  \n",
      "\n",
      "Instead we can use a `while` loop that continues to iteratively solve for the solution until the difference between two successive iterations is small enough.  \n",
      "\n",
      "How small is small enough?  That's a good question.  We'll try to work that out as we go along.  \n",
      "\n",
      "How can we best compare two successive potential fields?  The [L2 norm][1] is a good option.  It's defined as\n",
      "\n",
      "\\begin{equation}\n",
      "|\\textbf{x}| = \\sqrt{\\sum_{i=0}^n |x_i|^2}\n",
      "\\end{equation}\n",
      "\n",
      "There's one serious flaw with this.  We are effectively summing the difference at each point in our grid between successive iterations -- but what happens when the grid grows?  There are more points, in total, to compare and so more possible sources of error between the two solutions.  \n",
      "\n",
      "That doesn't seem fair, does it?  We can normalize the error by further dividing our existing L2 norm by the sum of the potential field.  \n",
      "\n",
      "For two successive iterations, the L2 norm is then calculated as\n",
      "\n",
      "\\begin{equation}\n",
      "|\\textbf{x}| = \\frac{\\sqrt{\\sum_{i=0, j=0}^n |p^{n+1}_{i,j} - p^n_{i,j}|^2}}{\\sqrt{\\sum_{i=0, j=0}^n |p^n_{i,j}|^2}}\n",
      "\\end{equation}\n",
      "\n",
      "The normalized `L2_error` function, then, will be\n",
      "\n",
      "[1]:http://en.wikipedia.org/wiki/Norm_%28mathematics%29#Euclidean_norm"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def L2_error(p, pn):\n",
      "    return numpy.sqrt(numpy.sum((p - pn)**2)/numpy.sum(pn**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's define a function that will do the iterations for us.  Three of the boundaries are Dirichlet boundaries and so we can simply leave them alone.  Only the Neumann boundary needs to be explicitly calculated at each iteration.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def laplace2d(p, y, dx, dy, l2_target):\n",
      "    '''Solves the diffusion equation with forward-time, centered scheme\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    p: 2D array of float\n",
      "        Initial potential distribution\n",
      "    y: array of float\n",
      "        Nodal coordinates in y\n",
      "    dx: float\n",
      "        Mesh size\n",
      "    dy: float\n",
      "        Mesh size\n",
      "    l2_target: float\n",
      "        Error target\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    p: 2D array of float\n",
      "        Potential distribution after relaxation\n",
      "    '''\n",
      "    \n",
      "    l2norm = 1\n",
      "    pn = numpy.empty_like(p)\n",
      "    iterations = 0\n",
      "    while l2norm > l2_target:\n",
      "        pn = p.copy()\n",
      "        p[1:-1,1:-1] = .25 * (pn[1:-1,2:] + pn[1:-1, :-2] + pn[2:, 1:-1] + pn[:-2, 1:-1])\n",
      "        \n",
      "        ##Neumann B.C. along x = L\n",
      "        p[1:-1, -1] = p[1:-1, -2]\n",
      "        l2norm = L2_error(p, pn)\n",
      "     \n",
      "    return p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Initial conditions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The problem states that our domain is set to zero everywhere as an initial \"guess\" except for \n",
      "\n",
      "$$p = \\sin \\left(  \\frac{\\frac{3}{2}\\pi x}{L} \\right) \\text{ at } y=H$$\n",
      "\n",
      "To define the domain, `numpy.zeros` will handle everything except that one Dirichlet condition."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "##variable declarations\n",
      "nx = 41\n",
      "ny = 41\n",
      "dx = 1.0/(nx-1)\n",
      "dy = 1.0/(ny-1)\n",
      "\n",
      "\n",
      "##initial conditions\n",
      "p = numpy.zeros((ny,nx)) ##create a XxY vector of 0's\n",
      "\n",
      "\n",
      "##plotting aids\n",
      "x = numpy.linspace(0,1,nx)\n",
      "y = numpy.linspace(0,1,ny)\n",
      "\n",
      "##Dirichlet boundary conditions\n",
      "p[-1,:] = numpy.sin(1.5*numpy.pi*x/x[-1])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's visualize the initial conditions using the `plot2D` function.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot2D(x, y, p)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The `p` array is equal to zero everywhere, except along the boundary $y = 1$.  Hopefully you can see how the relaxed solution and this initial condition are related.  \n",
      "\n",
      "Now run the iterative solver with a target L2 difference of $1\\tt{E}$$^-8$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = laplace2d(p, y, dx, dy,1e-8)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now try plotting this new value of `p` with the newly minted `plot2D` function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot2D(x,y,p)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Awesome!  That looks pretty good!  We'll need more than a simple visual check, though."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Relaxation & Convergence"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We're also going to write a function to \"reset\" the initial conditions for a given value of `nx` so we don't have to keep copying and pasting them.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def laplace_IC(nx):\n",
      "    '''Generates initial conditions for Laplace 2D problem for a \n",
      "    given number of grid points (nx) within the domain [0,1]x[0,1]\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    nx: int\n",
      "        number of grid points in x (and implicitly y) direction\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    p: 2D array of float\n",
      "        Pressure distribution after relaxation\n",
      "    x: array of float\n",
      "        linspace coordinates in x\n",
      "    y: array of float\n",
      "        linspace coordinates in y\n",
      "    dx: float\n",
      "        grid spacing in x\n",
      "    dy: float\n",
      "        grid spacing in y\n",
      "    '''\n",
      "    dx = 1.0/(nx-1)\n",
      "\n",
      "    ##initial conditions\n",
      "    p = numpy.zeros((nx,nx)) ##create a XxY vector of 0's\n",
      "\n",
      "    ##plotting aids\n",
      "    x = numpy.linspace(0,1,nx)\n",
      "    y = x\n",
      "    dy = dx\n",
      "\n",
      "    ##Dirichlet boundary conditions\n",
      "    p[:,0] = 0\n",
      "    p[0,:] = 0\n",
      "    p[-1,:] = numpy.sin(1.5*numpy.pi*x/x[-1])\n",
      "    \n",
      "    return p, x, y, dx, dy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Convergence Analysis"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to make sure that our Jacobi function is working properly.  Since we have an analytical solution, what better way than to do a convergence test?  As per usual, we will run our solver for several grid sizes and then check that the rate of decrease of error matches up with the increase in grid size.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx_values = [11, 21, 41, 81]\n",
      "l2_target = 1e-8\n",
      "\n",
      "error = numpy.empty_like(nx_values, dtype=numpy.float)\n",
      "\n",
      "\n",
      "for i, nx in enumerate(nx_values):\n",
      "    p, x, y, dx, dy = laplace_IC(nx)\n",
      "    \n",
      "    p = laplace2d(p, y, dx, dy, l2_target)\n",
      "    \n",
      "    p_an = p_analytical(x, y)\n",
      "    \n",
      "    error[i] = L2_error(p, p_an)\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.figure(figsize=(6,6))\n",
      "pyplot.grid(True)\n",
      "pyplot.xlabel(r'$n_x$', fontsize=18)\n",
      "pyplot.ylabel(r'$L_2$-norm of the error', fontsize=18)\n",
      "\n",
      "pyplot.loglog(nx_values, error, color='k', ls='--', lw=2, marker='o')\n",
      "pyplot.axis('equal');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That doesn't look like 2nd-order convergence.  What's going on?  The discretization is 2nd-order and Dirichlet conditions are order-agnostic (a  set value is a set value), but the Neumann boundary condition is 1st-order.  \n",
      "\n",
      "Remember when we said that the boundaries drive the problem?  One boundary that's 1st-order completely tanked our convergence.  Let's fix it!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "2nd-order Neumann Boundary Conditions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Up to this point, we have used the first-order approximation of a derivative to satisfy Neumann B.C.'s, e.g.\n",
      "\n",
      "\\begin{equation}\n",
      "\\frac{p^{n+1}_1 - p^{n+1}_0}{\\Delta x} = 0\n",
      "\\end{equation}\n",
      "\n",
      "which, solving for $p^{n+1}_0$ gives us\n",
      "\n",
      "\\begin{equation}\n",
      "p^{n+1}_0 = p^{n+1}_1\n",
      "\\end{equation}\n",
      "\n",
      "Using that Neumann condition will limit us to 1st-order convergence.  Instead, we can start with a 2nd-order approximation (the central-difference approximation):\n",
      "\n",
      "\\begin{equation}\n",
      "\\frac{p^{n+1}_1 - p^{n+1}_{-1}}{2 \\Delta x} = 0\n",
      "\\end{equation}\n",
      "\n",
      "That seems slightly problematic, since $p^{n+1}_{-1}$ doesn't exist.  But no matter -- according to the 2nd-order approximation,\n",
      "\n",
      "\\begin{equation}\n",
      "p^{n+1}_{-1} = p^{n+1}_1\n",
      "\\end{equation}\n",
      "\n",
      "Recall the finite-difference Jacobi equation with $i=0$\n",
      "\n",
      "\\begin{equation}\n",
      "p^{n+1}_{0,j} = \\frac{1}{4} \\left(p^{n}_{0,j-1} + p^n_{0,j+1} + p^{n}_{-1,j} + p^n_{1,j} \\right)\n",
      "\\end{equation}\n",
      "\n",
      "Notice that the equation relies on the trouble point $p^n_{-1}$, but according to Equation $(9)$, we have a value we can substitute, namely $p^n_{1}$, leaving us with the completed 2nd-order Neumann condition\n",
      "\n",
      "\\begin{equation}\n",
      "p^{n+1}_{0,j} = \\frac{1}{4} \\left(p^{n}_{0,j-1} + p^n_{0,j+1} + 2p^{n}_{1,j} \\right)\n",
      "\\end{equation}\n",
      "\n",
      "That's a bit more complicated than the first-order version, but it's relatively straightforward to code."
     ]
    },
    {
     "cell_type": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "$p^{n+1}_{-1}$ is *NOT* `p[-1]`.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`p[-1]` is a piece of Python code used to refer to the last element of a list or array named `p`.  $p^{n+1}_{-1}$ is a 'ghost' point that describes a position that lies outside the actual domain."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Convergence, Take 2"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here we can copy our old Jacobi function and then replace only the Neumann boundary condition line.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def laplace2d_neumann(p, y, dx, dy, l2_target):\n",
      "    '''Solves the diffusion equation with forward-time, centered scheme\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    p: 2D array of float\n",
      "        Initial potential distribution\n",
      "    y: array of float\n",
      "        Nodal coordinates in y\n",
      "    dx: float\n",
      "        Mesh size\n",
      "    dy: float\n",
      "        Mesh size\n",
      "    l2_target: float\n",
      "        Error target\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    p: 2D array of float\n",
      "        Potential distribution after relaxation\n",
      "    '''\n",
      "    \n",
      "    l2norm = 1\n",
      "    pn = numpy.empty_like(p)\n",
      "    iterations = 0\n",
      "    while l2norm > l2_target:\n",
      "        pn = p.copy()\n",
      "        p[1:-1,1:-1] = .25 * (pn[1:-1,2:] + pn[1:-1, :-2] + pn[2:, 1:-1] + pn[:-2, 1:-1])\n",
      "        \n",
      "        ##2nd-order Neumann B.C. along x = L\n",
      "        p[1:-1,-1] = .25 * (2*pn[1:-1,-2] + pn[2:, -1] + pn[:-2, -1])\n",
      "        \n",
      "        l2norm = L2_error(p, pn)\n",
      "     \n",
      "    return p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again, this is the exact same code, but now we're calling a Jacobi solver with a 2nd-order Neumann boundary condition.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx_values = [11, 21, 41, 81]\n",
      "l2_target = 1e-8\n",
      "\n",
      "error = numpy.empty_like(nx_values, dtype=numpy.float)\n",
      "\n",
      "\n",
      "for i, nx in enumerate(nx_values):\n",
      "    p, x, y, dx, dy = laplace_IC(nx)\n",
      "    \n",
      "    p = laplace2d_neumann(p, y, dx, dy, l2_target)\n",
      "    \n",
      "    p_an = p_analytical(x, y)\n",
      "    \n",
      "    error[i] = L2_error(p, p_an)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.figure(figsize=(6,6))\n",
      "pyplot.grid(True)\n",
      "pyplot.xlabel(r'$n_x$', fontsize=18)\n",
      "pyplot.ylabel(r'$L_2$-norm of the error', fontsize=18)\n",
      "\n",
      "pyplot.loglog(nx_values, error, color='k', ls='--', lw=2, marker='o')\n",
      "pyplot.axis('equal');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Nice!  That's much better.  It might not be *exactly* 2nd-order, but it's awfully close.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Learn More"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And for a detailed walk-through of the discretization of Laplace and Poisson equations (steps 9 and 10), watch **Video Lesson 12** on You Tube:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import YouTubeVideo\n",
      "YouTubeVideo('iwL8ashXhWU')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "References"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*  Thomas, James William. Numerical partial differential equations: conservation laws and elliptic equations. Vol. 3. Berlin: Springer, 1999.\n",
      "\n",
      "*  Thomas, J. W. (1999). Numerical partial differential equations: conservation laws and elliptic equations (Vol. 3). Berlin: Springer."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "css_file = '../../styles/numericalmoocstyle.css'\n",
      "HTML(open(css_file, \"r\").read())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}