{
 "metadata": {
  "name": "",
  "signature": "sha256:c334619ad18feb0b04faa57eab495e48215ad6b014d4ab8c09697baee5c35f56"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Multigrid Lecture and Homework Example"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Aron Ahmadia (US Army ERDC) and David Ketcheson (KAUST)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Teaching Numerical Methods with IPython Notebooks, SciPy 2014"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br /><span xmlns:dct=\"http://purl.org/dc/terms/\" property=\"dct:title\">This example</span> by <a xmlns:cc=\"http://creativecommons.org/ns#\" property=\"cc:attributionName\" rel=\"cc:attributionURL\">Aron Ahmadia and David Ketcheson</a> is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>.  All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "css_file = './example.css'\n",
      "HTML(open(css_file, \"r\").read())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
        "\n",
        "<style>\n",
        "\n",
        "@font-face {\n",
        "    font-family: \"Computer Modern\";\n",
        "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
        "}\n",
        "\n",
        "#notebook_panel { /* main background */\n",
        "    background: #888;\n",
        "    color: #f6f6f6;\n",
        "}\n",
        "\n",
        "div.cell { /* set cell width to about 80 chars */\n",
        "    width: 800px;\n",
        "}\n",
        "\n",
        "div #notebook { /* centre the content */\n",
        "    background: #fff; /* white background for content */\n",
        "    width: 1000px;\n",
        "    margin: auto;\n",
        "    padding-left: 1em;\n",
        "}\n",
        "\n",
        "#notebook li { /* More space between bullet points */\n",
        "margin-top:0.8em;\n",
        "}\n",
        "\n",
        "/* draw border around running cells */\n",
        "div.cell.border-box-sizing.code_cell.running { \n",
        "    border: 3px solid #111;\n",
        "}\n",
        "\n",
        "/* Put a solid color box around each cell and its output, visually linking them together */\n",
        "div.cell.code_cell {\n",
        "    background-color: rgba(171,165,131,1.0); \n",
        "    border-radius: 10px; /* rounded borders */\n",
        "    padding: 1em;\n",
        "    margin-top: 1em;\n",
        "}\n",
        "\n",
        "div.text_cell_render{\n",
        "    font-family: 'Arvo' sans-serif;\n",
        "    line-height: 130%;\n",
        "    font-size: 115%;\n",
        "    width:700px;\n",
        "    margin-left:auto;\n",
        "    margin-right:auto;\n",
        "}\n",
        "\n",
        "\n",
        "/* Formatting for header cells */\n",
        ".text_cell_render h1 {\n",
        "    font-family: 'Philosopher', sans-serif;\n",
        "    font-weight: 400;\n",
        "    font-size: 40pt;\n",
        "    line-height: 100%;\n",
        "    color: rgb(12,85,97);\n",
        "    margin-bottom: 0.1em;\n",
        "    margin-top: 0.1em;\n",
        "    display: block;\n",
        "}\t\n",
        ".text_cell_render h2 {\n",
        "    font-family: 'Philosopher', serif;\n",
        "    font-weight: 700;\n",
        "    font-size: 24pt;\n",
        "    line-height: 100%;\n",
        "    color: rgb(171,165,131);\n",
        "    margin-bottom: 0.1em;\n",
        "    margin-top: 0.1em;\n",
        "    display: block;\n",
        "}\t\n",
        "\n",
        ".text_cell_render h3 {\n",
        "    font-family: 'Philosopher', serif;\n",
        "    margin-top:12px;\n",
        "    margin-bottom: 3px;\n",
        "    font-style: italic;\n",
        "    color: rgb(95,92,72);\n",
        "}\n",
        "\n",
        ".text_cell_render h4 {\n",
        "    font-family: 'Philosopher', serif;\n",
        "}\n",
        "\n",
        ".text_cell_render h5 {\n",
        "    font-family: 'Alegreya Sans', sans-serif;\n",
        "    font-weight: 300;\n",
        "    font-size: 16pt;\n",
        "    color: grey;\n",
        "    font-style: italic;\n",
        "    margin-bottom: .1em;\n",
        "    margin-top: 0.1em;\n",
        "    display: block;\n",
        "}\n",
        "\n",
        ".text_cell_render h6 {\n",
        "    font-family: 'PT Mono', sans-serif;\n",
        "    font-weight: 300;\n",
        "    font-size: 10pt;\n",
        "    color: grey;\n",
        "    margin-bottom: 1px;\n",
        "    margin-top: 1px;\n",
        "}\n",
        "\n",
        ".CodeMirror{\n",
        "        font-family: \"PT Mono\";\n",
        "        font-size: 100%;\n",
        "}\n",
        "\n",
        "</style>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<IPython.core.display.HTML at 0x1074f8f10>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Multigrid\n",
      "\n",
      "Multigrid is one of the great algorithms for numerically solving PDEs.  It is also one of the few essentially optimal algorithms for solving a linear system of equations, since it computes the solution of an $N\\times N$ system in ${\\mathcal O}(N\\log(N))$ -- or even just ${\\mathcal O}(N)$ -- operations.  \n",
      "\n",
      "This notebook is meant to accompany a reading of [Section 4.6 of Randall LeVeque's text on finite difference methods](http://0-epubs.siam.org.library.kaust.edu.sa/doi/abs/10.1137/1.9780898717839.ch4).  Other good resources are cited there."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from matplotlib import animation\n",
      "from clawpack.visclaw.JSAnimation import IPython_display"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## A two-point boundary value problem\n",
      "Let's use Jacobi's method to solve the one-dimensional boundary value problem \n",
      "\n",
      "\\begin{align}\n",
      "u''(x) & = f(x) & 0<x<1\n",
      "\\end{align}\n",
      "\n",
      "with boundary conditions $u(0) = 1, u(1) = 3$ where\n",
      "\n",
      "\\begin{align}\n",
      "f(x) & = -20 + \\frac{1}{2} \\phi''(x) \\cos(\\phi(x)) - \\frac{1}{2} (\\phi'(x))^2\\sin(\\phi(x)) \\\\ \\\\\n",
      "\\phi(x) & = 20\\pi x^3.\n",
      "\\end{align}\n",
      "\n",
      "This problem appears on page 103 of the text. In the box below, we set up the grid, the right hand side,  and the boundary values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m=2**8-1                           # Number of grid points\n",
      "h=1./(m+1)                         # Mesh width\n",
      "x=np.linspace(0,1,m+2); x=x[1:-1]  # grid\n",
      "alpha=1.; beta=3.                  # boundary values\n",
      "phi = lambda x: 20.* np.pi * x**3\n",
      "f = lambda x: -20 + 0.5*120*np.pi*x * np.cos(phi(x)) - \\\n",
      "                0.5*(60*np.pi*x**2)**2 * np.sin(phi(x)) # RHS\n",
      "plt.plot(x,f(x)); plt.xlabel('x'); plt.ylabel('f(x)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Jacobi's Method\n",
      "A centered-difference method for this problem gives us the system of equations\n",
      "\n",
      "\\begin{align}\n",
      "\\frac{U_{i-1} - 2 U_i + U_{i+1}}{h^2} & = f(x_i).\n",
      "\\end{align}\n",
      "\n",
      "Solving for $U_i$ gives\n",
      "\\begin{align}\n",
      "U_i & = \\frac{1}{2} ( U_{i+1} + U_{i-1} ) - \\frac{h^2}{2} f(x_i).\n",
      "\\end{align}\n",
      "\n",
      "Jacobi's method consists of using the equation above as a fixed-point iteration:\n",
      "\n",
      "\\begin{align}\n",
      "U_i^{[k+1]} & = \\frac{1}{2} ( U_{i+1}^{[k]} + U_{i-1}^{[k]} ) - \\frac{h^2}{2} f(x_i).\n",
      "\\end{align}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the cell below, implement Jacobi's method for the BVP.  Don't forget about the boundary conditions.  Take 100 iterations with it and plot the result."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "U=np.linspace(alpha,beta,m)       # Just use a straight line as initial guess\n",
      "UU=[U]\n",
      "F=0.5*h**2*f(x)                   # Construct the RHS\n",
      "F[0]-=alpha/2.; F[-1]-=beta/2.    # boundary conditions\n",
      "e=np.ones(m-1)\n",
      "G=0.5*(np.diag(e,-1)+np.diag(e,1)) # Jacobi matrix (average neighbors)\n",
      "\n",
      "for i in range(100):\n",
      "    U=np.dot(G,U) - F\n",
      "    UU.append(U)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you can, try storing the solution after each Jacobi iteration in a list, and make an animated plot of the solution after each iteration."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xf = np.linspace(0,1,1000); \n",
      "uf = 1.+12.*xf-10.*xf**2 + 0.5*np.sin(phi(xf))  # Exact solution on fine grid\n",
      "u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))  # Exact solution on computational grid\n",
      "\n",
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, 1), ylim=(-3, 6))\n",
      "ax.plot(xf,uf)\n",
      "line1, = ax.plot([],[],'-')\n",
      "line2, = ax.plot([],[],'-o')\n",
      "plt.legend(['Exact','Jacobi','Error'],loc=2)\n",
      "\n",
      "def fplot(i):\n",
      "    line1.set_data(x,UU[i])\n",
      "    line2.set_data(x,UU[i]-u)\n",
      "    ax.set_title('# of iterations: '+str(i))\n",
      "    return line1, \n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=range(len(UU)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Convergence of Jacobi\n",
      "\n",
      "Notice how the high-frequency components of the error are damped very quickly, but the low-frequency error barely changes.  This is because the matrix $G$ appearing in Jacobi's method is tridiagonal, so in a single iteration each solution point only 'feels' the values of its immediate neighbors.  It takes very many iterations for low-frequency information to propagate.  As we'll see, part of the power of multigrid lies in enabling information to propagate rapidly between distant points in the solution.\n",
      "\n",
      "At each Jacobi iteration, the error in the $p$th Fourier mode  \n",
      "$$u^p = \\sin(p\\pi x)$$  \n",
      "is reduced by a factor of \n",
      "\n",
      "$$|\\hat{\\gamma}_p|= |\\cos(p\\pi h)|.$$\n",
      "\n",
      "Which modes of the error converge most slowly?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Under-relaxation\n",
      "The *under-relaxed Jacobi method* is obtained by taking, at each iteration, a weighted average of the previous iterate and the new iterate:\n",
      "\n",
      "\\begin{align}\n",
      "U_i^{[k+1]} & = (1-\\omega)U^{[k]}_i + \\omega\\left(\\frac{1}{2} ( U_{i+1}^{[k]} + U_{i-1}^{[k]} ) - \\frac{h^2}{2} f(x_i)\\right).\n",
      "\\end{align}\n",
      "\n",
      "Now at each iteration the $p$th Fourier mode is reduced by a factor of \n",
      "\n",
      "$$|\\hat{\\gamma}_p| := |1-\\omega + \\omega\\cos(p\\pi h)| = 1 + \\omega(\\cos(p \\pi h) - 1)|.$$\n",
      "  \n",
      "The code below plots this range of values. Try changing $\\omega$ to $\\frac{2}{3}$ and notice how the second half of the eigenvalues (for $p\\ge m/2$) are all shifted to have smaller absolute value.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M=30; H=1./(M+1)\n",
      "omegas=np.linspace(0.,1.)\n",
      "p=np.arange(1,M+1)\n",
      "\n",
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, M), ylim=(0,1))\n",
      "ax.set_xlabel('p', fontsize=20)\n",
      "ax.set_ylabel('$\\widehat{\\gamma}_p$', fontsize=20)\n",
      "line1, = ax.plot([],[],'o')\n",
      "\n",
      "\n",
      "def fplot(i):\n",
      "    omega = omegas[i]\n",
      "    gamma = abs(1.+omega*(np.cos(p*np.pi*H)-1))\n",
      "    line1.set_data(p,gamma)\n",
      "    ax.set_title('$\\omega = %02.2f $' % omega, fontsize=20)\n",
      "    return line1, \n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=range(50))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the cell below, implement under-relaxed Jacobi and run 100 iterations for the BVP with $\\omega=2/3$.  Does it make the results any better?  Why or why not?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "omega = 2./3.\n",
      "U=np.linspace(alpha,beta,m)  # Just use a straight line as initial guess\n",
      "UU=[U]\n",
      "F=0.5*h**2*f(x)\n",
      "F[0]-=alpha/2.; F[-1]-=beta/2.    # Construct the RHS, including boundary conditions\n",
      "e=np.ones(m-1)\n",
      "G=0.5*(np.diag(e,-1)+np.diag(e,1))\n",
      "niter=100\n",
      "\n",
      "for i in range(niter):\n",
      "    U=(1.-omega)*U+omega*(np.dot(G,U)-F)\n",
      "    UU.append(U)\n",
      "    \n",
      "    \n",
      "xf = np.linspace(0,1,1000); \n",
      "uf = 1.+12.*xf-10.*xf**2 + 0.5*np.sin(phi(xf))  # Exact solution on fine grid\n",
      "u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))  # Exact solution on computational grid\n",
      "\n",
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, 1), ylim=(-3, 6))\n",
      "ax.plot(xf,uf)\n",
      "line1, = ax.plot([],[],'-')\n",
      "line2, = ax.plot([],[],'-o')\n",
      "plt.legend(['Exact','Jacobi','Error'],loc=2)\n",
      "\n",
      "def fplot(i):\n",
      "    line1.set_data(x,UU[i])\n",
      "    line2.set_data(x,UU[i]-u)\n",
      "    ax.set_title('# of iterations: '+str(i))\n",
      "    return line1, \n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=range(niter))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Grid representations of a given frequency (*aliasing*)\n",
      "The code below plots a mode of a certain frequency, and then shows how it is represented on grids with decreasing numbers of points.  Notice how a *low frequency* mode (i.e., one with many points per wavelength) on a fine grid becomes a *high frequency* mode (one with very few points per wavelength).\n",
      "\n",
      "Eventually, on grids with very few points, the mode is *aliased* to a completely different mode."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import animation\n",
      "import matplotlib.pyplot as plt\n",
      "from clawpack.visclaw.JSAnimation import IPython_display\n",
      "\n",
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, 1), ylim=(-1.2, 1.2))\n",
      "\n",
      "p=5\n",
      "xf=np.linspace(0,1,1000) # fine grid\n",
      "\n",
      "line1, = ax.plot(xf, np.sin(p*np.pi*xf), '-r', lw=2)\n",
      "line2, = ax.plot([],[],'o-',lw=2)\n",
      "\n",
      "def fplot(m):\n",
      "    x=np.linspace(0,1,m+2); # grid\n",
      "    line2.set_data(x,np.sin(p*np.pi*x))\n",
      "    ax.set_title('m='+str(m))\n",
      "    return line2,\n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=(63,31,15,7,3,1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Restriction and Prolongation\n",
      "The multigrid idea is to solve successively on different grids, so that every mode has a small damping factor $\\hat{\\gamma}$ on some grid.\n",
      "\n",
      "In order to implement multigrid, we need a way to take a function with values on a fine grid and restrict it to its values on a coarser grid (restriction). We also need a way to take a function with values on a coarse grid and interpolate onto a fine grid (prolongation). The functions below do just this, using linear interpolation in the case of the prolongation operator. Make sure you understand what they are doing. You may find it helpful to look up slicing notation in the NumPy help online.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def coarsen(f):\n",
      "    return f[1::2]  # This slicing just takes the odd-numbered points\n",
      "\n",
      "def interpolate(f,alpha,beta):\n",
      "    m_coarse=len(f)\n",
      "    m_fine  =2*m_coarse+1\n",
      "    f_interp = np.zeros(m_fine)\n",
      "    f_interp[1::2]=f                          #Set the values of the odd numbered points\n",
      "    f_interp[2:-1:2]=0.5*(f[:-1]+f[1:])       #Set the values of the (interior) even numbered points\n",
      "    f_interp[0]=0.5*(f_interp[1]+alpha)     #Set the values of the endpoints\n",
      "    f_interp[-1]=0.5*(f_interp[-2]+beta)\n",
      "    return f_interp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Correction using a second grid\n",
      "\n",
      "Now we compute the residual $$r = f-AU$$ on the fine grid, coarsen it, and solve $$Ae=-r$$ for the approximate error on the coarse grid. Then interpolate and subtract the interpolated error from the fine grid solution.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#compute the residual\n",
      "A=2./h**2*(G-np.eye(m))\n",
      "F=f(x); F[0]-=alpha/h**2; F[-1]-=beta/h**2\n",
      "r=F-np.dot(A,U)\n",
      "print max(abs(r))\n",
      "omega=2./3.\n",
      "m2=(m-1)/2.               # Number of points in coarse grid\n",
      "h2=1./(m2+1)              # coarse mesh width\n",
      "x2=coarsen(x)             # coarse grid points\n",
      "r2=coarsen(r)             # residual restricted to coarse grid\n",
      "err2=np.zeros(m2)         # initial guess for the error\n",
      "e2=np.ones(m2-1); G2=0.5*(np.diag(e2,-1)+np.diag(e2,1))\n",
      "for i in range(3):                        #Solve Ae=-r by Jacobi iteration\n",
      "    err2=(1.-omega)*err2+omega*(np.dot(G2,err2)+0.5*h2**2*r2)\n",
      "err=interpolate(err2,0,0)     # interpolate the error\n",
      "print err[0]\n",
      "U=U-err                   # and use it to correct our solution\n",
      "plt.clf()\n",
      "plt.plot(x,u,x,U,x,U-u)\n",
      "plt.legend(['Exact','Jacobi','Error'],loc='best')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, in order to keep our multigrid code clean, we'll write a function to set up and take a fixed number of under-relaxed Jacobi iterations. The function returns the approximate solution and the residual.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Jacobi(U,f,alpha,beta,m,nu):\n",
      "    \"\"\"Perform nu Jacobi iterations on a grid with m points, with initial guess U, right hand side function f and \n",
      "       Dirichlet boundary conditions with values alpha and beta.  The function returns both the approximate\n",
      "       solution and the residual.\"\"\"\n",
      "    omega=2./3.\n",
      "    h=1./(m+1)\n",
      "    F=0.5*h**2*f.copy()\n",
      "    F[0]-=alpha/2.; F[-1]-=beta/2.\n",
      "    e=np.ones(m-1)\n",
      "    G=0.5*(np.diag(e,-1)+np.diag(e,1))\n",
      "    for i in range(nu):\n",
      "        U=(1.-omega)*U + omega*(np.dot(G,U)-F)\n",
      "    A=2./h**2*(G-np.eye(m))\n",
      "    FF=f.copy(); FF[0]-=alpha/h**2; FF[-1]-=beta/h**2\n",
      "    rr=FF-np.dot(A,U)\n",
      "    return U,rr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Multigrid V-cycle\n",
      "Now let's do a full V-cycle. Look carefully through the code below until you understand what it does.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "k=8;\n",
      "m=2**k-1\n",
      "rdep=k-1  # Recursion depth; this is how many grids down we want to go\n",
      "          # rdep=k-1 gives a full V-cycle\n",
      "nu=3      # Number of Jacobi iterations to take at each step\n",
      "U=np.linspace(alpha,beta,m)  # Initial guess\n",
      "x=np.linspace(0,1,m+2); x=x[1:-1]  # grid\n",
      "phi = lambda x: 20.* np.pi * x**3\n",
      "u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))\n",
      "f = lambda x: -20 + 0.5*120*np.pi*x * np.cos(phi(x)) - 0.5*(60*np.pi*x**2)**2 * np.sin(phi(x))\n",
      "alpha=1.; beta=3.\n",
      "F=f(x)\n",
      "r=[None]*(rdep+1); error=[None]*(rdep+1)        # This just initializes these lists to have the right length\n",
      "U,rr=Jacobi(U,F,alpha,beta,m,nu)  # Initial iteration on fine grid\n",
      "for i in range(1,rdep+1): # Going down\n",
      "    m=(m-1)/2 #  = 2**(k-i) - 1\n",
      "    r[i]=coarsen(rr)             # residual restricted to next coarser grid\n",
      "    error[i],rr=Jacobi(np.zeros(m),-r[i],0.,0.,m,nu)\n",
      "for i in range(1,rdep): # Coming up\n",
      "    m=2*m+1\n",
      "    err=error[rdep-i]-interpolate(error[rdep+1-i],0,0)     # Interpolate and subtract the correction\n",
      "    error[-i-1],rr=Jacobi(err,-r[rdep-i],0.,0.,m,nu)\n",
      "m=2*m+1\n",
      "U=U-interpolate(error[1],0,0)                   # final solution correction\n",
      "U,rr=Jacobi(U,F,alpha,beta,m,nu)                #Final iterations on original grid\n",
      "plt.clf()\n",
      "plt.plot(x,u,x,U,'o-',x,U-u)\n",
      "plt.legend(['Exact','Jacobi','Error'],loc='best')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since we have saved the errors on each grid, we can plot them to see how they represent the accumulation of errors at different scales."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, 1), ylim=(-5,5))\n",
      "line1, = ax.plot([],[],'-ok',lw=2)\n",
      "\n",
      "def fplot(i):\n",
      "    y = error[-i]\n",
      "    m = len(y)\n",
      "    x=np.linspace(0,1,m+2); x=x[1:-1]    \n",
      "    line1.set_data(x,y)\n",
      "    return line1, \n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=range(1,len(error)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Conclusion\n",
      "\n",
      "Now you have the basic ideas of multigrid.  Of course, for this 1-dimensional problem (where the system to be solved is tridiagonal) it does not provide a big advantage over more naive methods.  The power of multigrid is that it can be used in multiple dimensions, where the algebraic system is much more challenging for other methods.  Furthermore, multigrid can be applied to nonlinear and time-dependent problems."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Exercise #11:\n",
      "Modify the V-cycle code above to answer the following questions. Try to explain your results.  \n",
      "(a) How does the accuracy change as we change the number of Jacobi iterations performed at each step?  \n",
      "\n",
      "(b) Is it better to use a finer grid, or more Jacobi iterations if we want to improve the solution accuracy?  \n",
      "\n",
      "(c) What happens if we don't perform any Jacobi iterations in the \"up\" part of the V-cycle?  \n",
      "\n",
      "(d) What happens if we don't recurse all the way down to the 1-point grid?  \n",
      "\n",
      "(e) What happens if we use the original Jacobi method, or some other value of $\\omega$?  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Exercise #12\n",
      "Use the V-cycle code above to implement the full multigrid algorithm discussed in the text. You will probably find it helpful to write a Vcycle() function first.    "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "A Multigrid interact"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import display, HTML\n",
      "from IPython.html import widgets\n",
      "from IPython.html.widgets import interact, interactive\n",
      "\n",
      "def Jacobi(U,f,alpha,beta,m,nu,omega=2./3):\n",
      "    \"\"\"Perform nu Jacobi iterations on a grid with m points, with initial guess U, right hand side function f and \n",
      "       Dirichlet boundary conditions with values alpha and beta.  The function returns both the approximate\n",
      "       solution and the residual.\"\"\"\n",
      "    h=1./(m+1)\n",
      "    F=0.5*h**2*f.copy()\n",
      "    F[0]-=alpha/2.; F[-1]-=beta/2.\n",
      "    e=np.ones(m-1)\n",
      "    G=0.5*(np.diag(e,-1)+np.diag(e,1))\n",
      "    for i in range(nu):\n",
      "        U=(1.-omega)*U + omega*(np.dot(G,U)-F)\n",
      "    A=2./h**2*(G-np.eye(m))\n",
      "    FF=f.copy(); FF[0]-=alpha/h**2; FF[-1]-=beta/h**2\n",
      "    rr=FF-np.dot(A,U)\n",
      "    return U,rr\n",
      "\n",
      "def Vcycle(k=8,nu=3,omega=2./3,smooth_after_interp=True):\n",
      "    m=2**k-1\n",
      "    rdep=k-1  # Recursion depth; this is how many grids down we want to go\n",
      "              # rdep=k-1 gives a full V-cycle\n",
      "    alpha=1.; beta=3.\n",
      "    U=np.linspace(alpha,beta,m)  # Initial guess\n",
      "    x=np.linspace(0,1,m+2); x=x[1:-1]  # grid\n",
      "    u = 1.+12.*x-10.*x**2 + 0.5*np.sin(phi(x))\n",
      "\n",
      "    F=f(x)\n",
      "    r=[None]*(rdep+1); error=[None]*(rdep+1) # This just initializes these lists to have the right length\n",
      "    U,rr=Jacobi(U,F,alpha,beta,m,nu,omega)  # Initial iteration on fine grid\n",
      "    for i in range(1,rdep+1): # Going down\n",
      "        m=(m-1)/2 #  = 2**(k-i) - 1\n",
      "        r[i]=coarsen(rr)             # residual restricted to next coarser grid\n",
      "        error[i],rr=Jacobi(np.zeros(m),-r[i],0.,0.,m,nu,omega)\n",
      "    for i in range(1,rdep): # Coming up\n",
      "        m=2*m+1\n",
      "        err=error[rdep-i]-interpolate(error[rdep+1-i],0,0)     # Interpolate and subtract the correction\n",
      "        if smooth_after_interp:\n",
      "            error[-i-1],rr=Jacobi(err,-r[rdep-i],0.,0.,m,nu,omega)\n",
      "        else:\n",
      "            error[-i-1] = err\n",
      "\n",
      "    m=2*m+1\n",
      "    U=U-interpolate(error[1],0,0)                   # final solution correction\n",
      "    U,rr=Jacobi(U,F,alpha,beta,m,nu,omega)                #Final iterations on original grid\n",
      "    \n",
      "    fig, ax = plt.subplots(figsize=(8,6),subplot_kw={'axisbg':'#EEEEEE',\n",
      "                                       'axisbelow':True})\n",
      "    ax = plt.axes(xlim=(0, 1), ylim=(-5,6))\n",
      "    ax.plot(x,u,x,U,'o-',x,U-u)\n",
      "    plt.legend(['Exact','Jacobi','Error'],loc='best')\n",
      "\n",
      "interact(Vcycle, k=(4,10,1), nu=(1,4,1), omega=(0,2,1./6))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}