PS3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast PDE/IE course, Skoltech, Spring 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Set 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral Methods\n",
    "\n",
    "Consider a problem of solving\n",
    "\\\\[\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f, \\\\\n",
    "u|_\\Gamma &= 0.\n",
    "\\end{align*}\n",
    "\\\\]\n",
    "by a spectral method, where $\\Omega = [0,1]^2$ is a unit square and $\\Gamma=\\partial\\Omega$.\n",
    "\n",
    "In a spectral method you will be working with a representation of the solution in the form\n",
    "$$\n",
    "u(x,y) = \\sum_{i=1}^N \\sum_{j=1}^N \\hat{u}_{i,j} \\sin(\\pi i x) \\sin(\\pi j y),\n",
    "$$\n",
    "where $\\hat{u}^N$ is just an $N\\times N$ array.\n",
    "\n",
    "Knowing $u(x,y)$, it is often practically not possible to compute $\\hat{u}_{i,j}$ exactly. Hence one needs to use the [Discrete Sine Transform](http://docs.scipy.org/doc/scipy-dev/reference/tutorial/fftpack.html#discrete-sine-transforms) (DST).\n",
    "The python's <tt>dst</tt> computes a one-dimensional DST, but we need a two-dimensional (2D) DST, which we compute in the following way:\n",
    "    \n",
    "* Compute values $u(x,y)$ for $N^2$ points inside the domain $(x,y) = (h k, h \\ell)$, $1\\leq k,\\ell \\leq N$, $h = 1/(N+1)$.\n",
    "\n",
    "* Apply <tt>dst</tt> to each line in this array. This way you will be computing a \"partial\" 2D DST\n",
    "$$\n",
    "u(x,y) = \\sum_{j=1}^N \\tilde{u}_{j}(x) \\sin(\\pi j y),\n",
    "$$\n",
    "\n",
    "* Apply <tt>dst</tt> to each column in the resulting array. This way you will be computing a full 2D DST.\n",
    "<br><br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1 (Spectral Methods) (60pt)\n",
    "\n",
    "* Part (a)\n",
    "    - Take $u^0(x,y) = \\sin(\\pi x) \\sin(\\pi y)$ and compute $f = -\\Delta u^0$. \n",
    "    - **(6pt)** Calculate, explicitly, the coefficients $\\hat{f}_{i,j}$ in \n",
    "    $$\n",
    "    f(x,y) = \\sum_{i=1}^N \\sum_{j=1}^N \\hat{f}_{i,j} \\sin(\\pi i x) \\sin(\\pi j y).\n",
    "    $$\n",
    "    Then, using values of $f$ at $N^2$ points, calculate the coefficients $\\hat{f}^{N}_{i,j}$ by the procedure outlined above.\n",
    "    (Here $\\hat{f}$ is an \"infinite array\" that you will compute explicitly, while $\\hat{f}^N$ will be an NxN array that you will compute using python.)\n",
    "    Compare $\\hat{f}$ and $\\hat{f}^N$.\n",
    "    - In this case, obviously, a solution to the problem is $u=u_0$, but let us now pretend we do not know it.\n",
    "\n",
    "    - **(6pt)** Compute the coefficients $\\hat{u}^N$ from $\\hat{f}^N$ found in part (a).\n",
    "\n",
    "    - **(6pt)** We will estimate the error between the two solutions in the following way. We will take $N^2$ points in our domain, of the form $(h k, h \\ell)$, same as before. We will then be interested in\n",
    "    $$\n",
    "    {\\rm err}_N = \\max_{1\\leq k,\\ell\\leq N} |u^N(h k, h \\ell) - u^0(h k, h \\ell)|\n",
    "    $$\n",
    "    which we call an error of the solution. Calculate the error of your solution for a number of values of $N$. Explain your results\n",
    "\n",
    "* Part (b)\n",
    "    - Let us now take something a little more complicated,\n",
    "    $$\n",
    "    u^0(x,y) = \\sin(\\pi x^2) \\sin(\\pi y^2),\n",
    "    \\qquad\\text{and}\\qquad\n",
    "    f = -\\Delta u_0,\n",
    "    $$\n",
    "    - You cannot compute the coefficients $\\hat{f}$ explicitly, so you'll have to live with only $\\hat{f}^N$.\n",
    "    - **(6pt)** Compute $\\hat{u}^N$, the error ${\\rm err}_N$, and report the error for different values of $N$. Would you say the error decays fast as $N$ increases?\n",
    "    - Suppose that the spectral method has an order of convergence, in other words the error behaves like ${\\rm err}_N = C N^{-\\rm ord}$ and you need to find ${\\rm ord}$.\n",
    "        -  **(6pt)** Derive the formula\n",
    "        $$\n",
    "        {\\rm ord}_N = \\frac{\\ln({\\rm err}_N)-\\ln({\\rm err}_{2N})}{\\ln(2)}\n",
    "        $$\n",
    "        -  **(6pt)** Hence compute ${\\rm ord}_N$ for $N=1,\\ldots,20$ and comment on your results. (You should start with taking each value of $N$ between 1 and 20 to understand the behavior, but you don't need to present all these numbers in your report, as long as you can illustrate the right behavior.)\n",
    "\n",
    "* Part (c)\n",
    "    - Let us now play a \"fair game\": take\n",
    "    $$\n",
    "    f(x,y) = 1.\n",
    "    $$\n",
    "    We then do not know the solution.\n",
    "    - Use your code from part (b) to compute $\\hat{f}^N$, and $\\hat{u}^N$.\n",
    "    - **(6t)** Instead of the exact error we have to use the error estimate\n",
    "    $$\n",
    "    {\\rm errest}_N = \\max_{1\\leq k,\\ell\\leq N} |u^N(h k, h \\ell) - u^{2N+1}(h k, h \\ell)|\n",
    "    $$\n",
    "    (Why did we take $2N+1$ instead $2N$?)\n",
    "    Report the values of ${\\rm errest}_N$ for a sequence of values of $N$.\n",
    "    \n",
    "    -  **(6pt)** Hence compute ${\\rm errest}_N$ for a sequence of $N$ and comment on your results.\n",
    "    -  **(6pt)** Finally, using the same formula (but with ${\\rm errest}$),\n",
    "    $$\n",
    "    {\\rm ord}_N = \\frac{\\ln({\\rm errest}_N)-\\ln({\\rm errest}_{2N})}{\\ln(2)}\n",
    "    $$\n",
    "    compute ${\\rm ord}_N$ for a squence of values of $N$. Comment on your results. \n",
    "*  **(6pt)** Compare the behavior for ${\\rm err}$, ${\\rm errest}$, and ${\\rm ord}$ for parts (b) and (c). What is the main reason for the qualitative difference in the speed of convergence?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2 (Multigrid) (40 pts)\n",
    "Consider a Poisson equation\n",
    "\\\\[\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f, \\\\\n",
    "u|_\\Gamma &= 0,\n",
    "\\end{align*}\n",
    "\\\\]\n",
    "where $\\Omega = [0,1]^d$, $d = 2, 3$ and $\\Gamma=\\partial\\Omega$.\n",
    "\n",
    "Although the multigrid method has optimal complexity, there are some precomputations to be done. Therefore, it may have a bigger constant than, for instance, the spectral method or even sparse LU. In this problem you will be asked to find out in which case (2D, 3D, grid sizes) multigrid method is more appropriate.\n",
    "\n",
    "* (**20 pts**) Implement V- and W-cycles of multigrid method for $d=2$ and $d=3$. Assume that the equations are discretized with finite differences with the standard second order \"cross\" stencil.\n",
    "**Note:** all operations should be implemented with linear complexity. Take a smoother of your choice.\n",
    "\n",
    "\n",
    "* (**5 pts**) Set $N=128$. Choose $f$ such that you know analytical solution. Plot errors as a function of the cycle number for the implemented V-,W-cycles of GMG and for AMG (from PyAMG), $d=2$ and $d=3$. Which approach requires fewer number of iterations? Which approach is faster? **Note:** when doing plottings do not forget to inculde title and axes names\n",
    "\n",
    "<!---\n",
    "* How many iteration for your smoother are required to reach the same accuracy as for multigrid?\n",
    "-->\n",
    "\n",
    "* (**5 pts**) Implement a spectral method for $d=3$. Implement fast Poisson solver for $d=2$ and $d=3$: the only difference with the spectral method will be eigenvalues \n",
    "$$\\lambda_{mn} = -\\frac{4}{h^2} \\left( \\sin^2 \\frac{\\pi m h}{2} + \\sin^2 \\frac{\\pi n h}{2} \\right),$$\n",
    "(eigenvalues of the discrete problem ($d=2$) with the \"cross\" stencil) instead of \n",
    "$$\n",
    "\\lambda_{mnl} = -\\pi^2(m^2 + n^2)\n",
    "$$\n",
    "(exact eigenvalues of the initial continuous problem).\n",
    "\n",
    "\n",
    "* (**10 pts**) Set $f\\equiv1$. Investigate when (for what cases) GMG is more appropriate than the sparse LU, the spectral method or the fast Poisson solver. To do so compare timings for different $N$, $d$ and accuracies. Make conclusions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "breakpoint": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<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=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 100%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 100%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 100%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/alex.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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
}