Lecture-19-Spectral-Solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 19:  Numerical Solutions to the Diffusion Equation  (Spectral Method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## What to Learn?\n",
    "\n",
    "* The assumptions of the spectral method\n",
    "* The possible choices for basis vectors\n",
    "* The linear independence of resulting ODEs\n",
    "* One way to code the FFT/IFFT\n",
    "* Finite differencing of the amplitudes\n",
    "* Useful data structures and pointer swapping using tuples\n",
    "\n",
    "## What to do?\n",
    "\n",
    "* Solve the diffusion equation using a spectral method.\n",
    "* Analyze two example implementations and determine the difference scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Introduction\n",
    "----\n",
    "\n",
    "In this course the term _spectral method_ refers to a numerical solution composed of a finite set of basis functions and the time evolution of their amplitudes.  Practically, the choice of basis functions will satisfy the boundary conditions and the initial amplitudes will satisfy the initial conditions.  The evolution of the amplitudes can then be written as a finite difference in time.  Although other choices for basis functions are [possible](./images/Boyd_Table.png)$^*$, we will focus on Fourier functions.\n",
    "\n",
    "$*$ Chebyshev & Fourier Spectral Methods, Boyd, John P., Springer 1989"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Developing the Logic of the Spectral Method\n",
    "----\n",
    "\n",
    "The spectral method assumes that the solution to a PDE (the function $c(x,t)$) can be represented as a series expansion that contains time dependent amplitudes (the $a_k(t)$ factors) and spatially varying functions (the $\\phi_k(x)$ terms are unit Fourier vectors).  Note the similarity to the approach taken in the separation of variables method in the following discussion.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To begin, we write:\n",
    "\n",
    "$$\n",
    "c(x,t) = \\sum_{k=0}^N a_k(t)\\phi_k(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Choosing the unit vectors from an orthonormal set (like Fourier) permits decoupling the summation into $N$ independent equations that can be solved simultaneously.  Furthermore, the assumption of the orthonormal set allows the time dependence to be placed within the amplitudes so that the independent variables are seperated.  Therefore we refer to the $a_k(t)$ as the set of amplitudes with one $a_k$ for each of the $k$-basis vectors, $\\phi_k(x)$.  There is no restriction on the value of $N$ although practical considerations that relate to diffusive problems and the \"smoothness\" of the solutions will require just a few $N$ terms to achieve acceptable accuracy.  The form of $\\phi$ depends on the problem being solved; the boundary conditions and the initial conditions will affect this choice.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our example problem:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c(x,t)}{\\partial t} = \\frac{\\partial^2 c(x,t)}{\\partial x^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "with the initial and boundary conditions:\n",
    "\n",
    "$$\n",
    "c(0,t) = 0\\\\\n",
    "c(L,t) = 0\\\\\n",
    "c(x,0) = c_0(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The problem requires that our boundary values for the $c(x,t)$ at $0$ and $L$ are zero.  One possible choice is a series of $\\sin$ terms.  This satisfies the boundary conditions and enforces perodicity of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "\\phi_k(x) = \\sin \\left( \\frac{k\\pi x}{L} \\right)\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Begin by importing relevant libraries and defining the symbols we expect to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "x, y, z, t = sp.symbols('x y z t')\n",
    "k, m, n = sp.symbols('k m n', integer=True)\n",
    "f, g, h = sp.symbols('f g h', cls=sp.Function)\n",
    "sp.var('a_k, phi, c', cls=sp.Function);\n",
    "sp.var('L', real=True);\n",
    "sp.init_printing();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The assumed form of the solution is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOAAAAAmBAMAAAAvsop7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXYyIt1Uu83vq5kQZkTr643dAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEZklEQVRYCbVWXYgbVRT+JjP5nyRrVSgo7GypUHSLwYJoUXdeyoJYN2JFHwTDqoggbBBBKyyE+lBR1GiLLyoM22JLVVwFofTF0VKq+JB1QahIMQqCPm2isl2728Zz7tw7c2dMfJoeyJzvfPeb78xM7pwEuBZRWRnlanuj2FS4XeRy+3A43IS94/iv25Tn7wqknavvA8b18/t3uvgKPziPKf+Mp1DK+ds6cCduyFC6GZP+U6H9GyFKF5wSdq2SD/g4rHnPEXMNwvybTU0vK7z/1DqUPK1IDxaEbXml5pOn+ZdmbG5pRXpwzmWvrl/rAc9l+pjmMohDCqSavxBut6A4gep6vm9pL+WPYaP9IQpAxU0QiVLox4jeEdpFlCZgL03fdlE7da4nC6ulsQIeTxKxWupDkRk0ERpL/9bkWfbBhkB000EUXAmCdM7FNzEiWbBeF5mDSFHpRzhEhWUBS23J/BSuCLDgw3AEGnNg/ThRObiZ+Jm1pqhpCwVxQGaZTgN2J07FK9aPE2Vbca2o1nyRcuIVBazERT1Eq28KxeiD0Osia08krC1HOEQ3Bshcp7xtRz1DmmiuV09cWgJuCsUSrE45iiJ9JDJ3r7ZnL+Ou9+64VQhqoUzJKb+O+WOU7CvAvma+nV+BPtdzfVq7wHLjVRF8zbkDOIP7PWbB+lB0HmjgNeBRxxIPrEuLybDXMw9SL4AG3dsoNbPN2Fwvt2npafrokTuCNT/XERTrQ9EfPhy8DBwEhOfciIa5rWlb0AOUebwVe7G5np0gboY+sTCnvnTNoCHrQ1H20k5ww5fE9dOJTT6Nfn7D6KO8yRzFJootSmwg5Qy6PTokG9qH3Zm61jAU2duHy9zwRdVQ3Ar7RJFdej4oNtHlu+FHpM31NdjAb6wwHhHhEcw3MHPP2U7+MmHWh6IHULmqNxz1SLu9YKjZA9Toepq8CbS5vh1Pyk1D5ipqE1i4b0/nvEsE60PRPPCJ3rDrqFOi/DO6Tb4z6wrKHr6rG7TNg7nebRN9xnaAzwjoQc9+8t6zR9vMkT4SPevjBX6e6pHSlf0nTqBWZ5rfw+suOqg21FzPPk70LA/4VxJnWR9NfX3IHLzLNOkj0S+r3/f2Dj/eOzz2zD8eLRb5kIhFZD73iatelQs8qoJQ/3iSwydYNRt3u4ykfqQo32YFR+XTjQ8DpI6VLYl4GAehnofhSCKWzIbRZkLqDYeLRBidkKjRk4hFuJZ3Jc/fnojRP0+zg8yGS+tSP1KUU3dBO1ldfuBJ29uTyGpJUJGZ9sH/hdSPFPG8lDHZU0jmaOyJvwzaasXVihFQ6MeIaAjIeEsBlRdchdLM0T+ljaTt6SSRSt1tSpuq9scjoD5IpUHSJNOWjKHeOqXIRBtYUankI9Kl0Cfga5bFFa1IEV6QXjXa+7GNtSvFJrpVwQmqNcplbcE6qhVpQvvhwG2yDuzWjM81tSJVuM8Vdk/Q70o0doDFVJvoZhbdGuyTw1Mnh17E2z3C/wIRcCmvaLYQXQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$c{\\left (x,t \\right )} = \\operatorname{a_{k}}{\\left (t \\right )} \\sin{\\left (\\frac{\\pi x}{L} k \\right )}$$"
      ],
      "text/plain": [
       "                    ⎛π⋅k⋅x⎞\n",
       "c(x, t) = a_k(t)⋅sin⎜─────⎟\n",
       "                    ⎝  L  ⎠"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "elementK = sp.Eq(c(x,t),a_k(t)*sp.sin(k*sp.pi*x/L))\n",
    "elementK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOBAMAAADkjZCYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXYyEM1Embsi72ZU3au6f2Q3AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAVUlEQVQIHWNgYBBUNGBgcE1gD2BgaGfgaGBg+8kABBwLQCR/sFAZA4N/NAPXAQb/LwzsSxj4LzBwf2PgDWDg/s3AtwDEZgaKAxU/ZPA6wMDAWP6YAQCADRI3fgJiQgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$6$$"
      ],
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.acos()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Proceed by substituting the series expansion into the PDE and performing the differentiations as defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANEAAAAvBAMAAACPjDK1AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMnaZiVQiZrurRO9IrJcIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE/klEQVRYCbVXXWgcVRQ+szOb2dmf7FgNgkZ3H1QUiglI0WDFwVLIg+yu0VSiwax/L/Whq2ikgu2CCAWhWbQv1p+MDwotQqJSoUHYrSC+FLpgwRcxKQriHya4pq2mrOfcn7mzM5N0UpIDO/c73zn3uzNn7t6zCxAw84WnBaNQIGWL3IPwtVBSaIukAzIjUG5wSqFAyha5D8Bsm0sptEXSIZkZW1IKSWZrx10AN//EJBFtp5lVVL+DVmCIwDbZJ6R7iS4MEdgey9VyDuRWUJyh7VmEqT505C0bzGHEDG1ypZOR+YMBdqrb7a7CX90uQKKamnQYCuRcxU0WAZIo063CwOQHj1dEeg5pn5kvv3l8osGJQn1i0fXF4sIf8Qvy6ujEx++CVkvv0oty3qQEbByFg0lHMDMD9k09wXiO9jOAUTcbj+G+tc2VlLhvrJHbKzCfxVtidqAqwOaGbIvyM3AGQIPEvJps0Q5TZrR06ex+siLhZkZ+Yj4MwzQp75c42yPT5+ZtTmhXrCs9oZjOTpa3H9gTlB3frLJQ5lTBzrc5Si9rl0xfXkxoXaZEbQXwNlP2FCSVfKLi13gU+mvcz85D5xpW6hum6clh+AVgsd6E40o9vaQwwAJkxUoZF/Y94g8pbLmETyjCh/rZjePkbwEGxwf2uyp2DW/jPTY725AiJabOvVlHsmoUXYEec5N2F76IEdCKclpZlIH8OUn6R9YVwPvJ4A9tiNNVfNm4tV6Jyno2imRdAWAqKrYRp2PZ+ooAt0Yl3R9Bsq6A/IxX74ikKOodJHUsWKItoscclbaqoIdYV0CvXPeo9UF6XNVqH8APza8qvB3vPTlmNVvG3c+Nf8Nm/xOhwboC8gU3Ihik3mj0tSxxnCxgcAY/6SLujHnQIdEC/RTMUm20f/ESNNkVes6mYJL0jzpY7KPcuw+HZ/CTGsaefAGSkGlBZgnyLnJW1EqyK+RrsOMestsAqIGFbBkVYO/0ZSgRAPgTP/fix6DAg6cabKUiP00Nts3yfg2AA1VMRMtX2bDxZcdpo+NbSaPDl61049AqW2lerBT1TLIrxFrpCzA615Wg2cIVsHq5ZbAhh9WzbJhrU/XESlHV87pCobbO40x/aXuRi5DqnC9pp4nAHWEWUw7bEakKZF3fSlrE3vO6wqzr6fUAs5EtesTvYHY+Kr3N/Gn84rYG8IK3mLoAmQY+kHwmuOhN8YDXFdb7PukNg/UaNuP9518cq5QP1cn5Dou24OI318EyvvbE68a5tRvOrY02/wA0diGgzOsKi21FCvRpd2cd9HbwmN9j/EoJdBqR0VkRss89hol4HoE59TY8/m+G6Bj12x54iVw6YcnUqcF9dj2kMBdRPvhiks3xausVSbDRatamsGEAsH0BWhFhyGa8+xYivozdPiygyX+r4KuPNN4JM6KIvSkFV/pCRLrRJ1VimeKySipXIMslcEJ4vQPueGFcRHo4ptQG89h8keD3cN5jYoP0kkzlItLDkXpa0BZryBh3HvksGLi6r4n+ir+MSMRv/a7f4/hwG8f+bpcVMRzekBmTUSYiHRqf8jsCD+GoRfBxqH5HZIVEfouY/h9y0S88IjlAUX9kFhSh/1VBy9ExWQ2ycf3beWJI5MNGWMHs4O9dN8zHY653WF5QRKNWHTQdd8JhO8jG9a0RlhkUwXM2ZINDa2eH1kJ0bOIYZYZEbiH2fzHBRasZ8zVgAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- \\frac{\\pi^{2} k^{2}}{L^{2}} \\operatorname{a_{k}}{\\left (t \\right )} \\sin{\\left (\\frac{\\pi x}{L} k \\right )}$$"
      ],
      "text/plain": [
       "  2  2           ⎛π⋅k⋅x⎞ \n",
       "-π ⋅k ⋅a_k(t)⋅sin⎜─────⎟ \n",
       "                 ⎝  L  ⎠ \n",
       "─────────────────────────\n",
       "             2           \n",
       "            L            "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spaceDeriv = elementK.rhs.diff(x,2)\n",
    "spaceDeriv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKYAAAArBAMAAAAAiO5CAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMnZmzRC73UTvIomZVKu7zOipAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEZklEQVRIDZ1WXYgbVRT+JsnkZ5Jstv6hgs1ouyhd0bBdKKJgXKH4UHQqpaUIJrhF2D5oBPGhggYRpPiwqUUflgoDxdUH6aalIoiwQQo+FNq4oKIvO421D0Vr1N1S3abx3HNnMneSmM3uIdz7ne8753Bncs9JgE3ZXSc2lTY4qTVY3owaX91M1uCcZGmw3qOGcz0UEbqjsJGa4gwDt1PQo+12+yb0B66+c8ZLedcD+PxMtt5xhgLxD4Hk6aWTO2y8gF3m615SyHFRfA6vWh473H5wFLgPX4Vo+wKL1s+drKdddDCHcx1yOHCNw4ojdBQLu5WcBfdwFy28qdBDQO1fEaQ5MY79S8kYcaTzEfQbCj0EjHBmKpcRp9L+VDI0eSn1NaSbCj0EXLBFUMHKlIHvQ02MCVfaFG96C9G8/5Y9cdD+LYsPIlFFvBVtGjk/+DEJp7DLqfrsEOgDjpnGSBV6Y+zhnUrKQpmdF8fvbUikaARPCjdsi5XemqzD2FDfIDPUQjMlRnT0AWYUWbwqQ7Q1PzTc9HEHua04ku8wfUDEBg7Y2NcrpeSRgkKmzj59YwPsCdKWLSTNnphYsYcCZi0m03x1++hMHaGV+kuvsGds8eMyNR930JcSaYMmsVGiINFfz0DbNpGfv4U7Ttz2iCkyM7wKpNgeLF0iV/9b4RhOjJseFaohfuV6A/ga+4ESngXeMA1+skLOC/J3vRX6jMupnSrk9BFcwCmHI6OUyP01iaMWTLwCzACcttCnZnp1TGdauR9cJ32MXnW6wjhWB1J5gj8idn0HRM1DAJ+i37RN3eQswNtdl272+EVbkzUTZSBWJSUL/Wy7Jmq+7NXMiQSa8h1rItZ4SZCypsJD321nR5WahTIFZfEJwjfUmv2evVCWXUkTKWjRErJ37q1EbxEtnn0WOvAeloBf1JoFM5gmvCdRqIunMrq/90wVy3dvqey3SRTf0Vl8B0ziJwuHxYN7z05hPXYFmVFB99zPRBGLt++dy4uMZA24oJvAebw/8Xh5a/vXre1LP/zjkJQQS5dNI/SNRVy8e7wbv40/P6WtfSzi4yVgfieBp4QXsGg+4AaccN9/H1rpHluEid4k436S0FuTFQ+F/1j53cNy9zRtpqYIWimZF66YIWR9ZkjaP0uGniZgMcd1aU7RVHNtfi20IpyoWNBv1ik9PVvlIH/p9C3NKZpqXWYUmTjcRQuX7r9ri2UPufuy7QKaU+e6NHIDvx0B2f0hI+65AE9Op84eOdW69f/1C3VPWvGAtx9nYJy/vyWnmsevu4fybki8uxNDFVaO2tHVjf5rOObWTHbf8EROKPpxMdLosxGbdIMjTQKWkrmdcWpVjDSeaoq2DoyYMiBT9P8BCMaYY556nEYafTZi+msyetakwa0kHqizU6BZ5E41RVwPfmpzxOIosE2JnZaY7sWhtM1TTRHXgwYVI3ubetdvVOhlZhExw28l5VSTxPCrfrl97XLb6U3QG6dPPSSnWq/IzH+AkiJJ7b8wLgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\sin{\\left (\\frac{\\pi x}{L} k \\right )} \\frac{d}{d t} \\operatorname{a_{k}}{\\left (t \\right )}$$"
      ],
      "text/plain": [
       "   ⎛π⋅k⋅x⎞ d         \n",
       "sin⎜─────⎟⋅──(a_k(t))\n",
       "   ⎝  L  ⎠ dt        "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timeDeriv = elementK.rhs.diff(t,1)\n",
    "timeDeriv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Our final differential equation represented in $a(t)$ is therefore:\n",
    "\n",
    "$$\n",
    "\\sum_{k=0}^N \\sin{\\left (\\frac{\\pi x}{L} k \\right )} \\frac{d a_k{\\left (t \\right )}}{d t}  = - \\sum_{k=0}^N \\frac{\\pi^{2} k^{2}}{L^{2}} a_k{\\left (t \\right )} \\sin{\\left (\\frac{\\pi x}{L} k \\right )}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a reminder, a $\\sin$ series is orthogonal over $0 < x < 2\\pi$ if the following integral is zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMQAAABKBAMAAAAWMubwAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiZm7ECKrzd0yRFR2Zu/H3ntdAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD70lEQVRYCe2YXWgcVRTH/7Mzk52P3e4YpFSpNAVBfQgsRREkdm9fLELpTikVP5AsSPElkqWlDyrYCH482JJ9lYq7vrVQmn2oIAjNpA0VSyR5Ux9KUmh9sBVX3UZMrNszs5mZ3cmdPO0pfXAgueee+Z3zvx9z585dABjeI+g/56W4Q1XO/JTbcM273BINfY1ZArDb+slO5x9OHcPLfD7/pOCUGIbIuSOcCnnKruqsEl/hHA5bFcZemO9/+yxGzRVGiW2dTgsf4yNGif9TP2QjkHlEcLfo6sEKs4S2ZnjMErkWswBtS0V2iXHuYQIWG+y9KAl2if/YFcxVdoncv+wSKv9XTqHN3ov6FLvEckUqoS8l3Hk34Qiq/ZicwaQXhX4ZmzgceUPjaGhsgYWMuXckpKmcbkSVuhOZOBGbG9b+0LMFFjGFZkgDWkfEldjyP+ESl+IlHFRNYhEzI2LYli9uY/PIaytxVGglsYjZERJUKn/Flddqkf1qZMXGexvmVljIHIvDkP2jW/ni+nd6qWJ+8MzChO+4DFjbjyxMWDtfAA5cfX3Oo0Mb+VOwPgb67u/XoS0cPb/TT4XCelBoRRjIVmCMYdwfoxepf+oYytdwCZqTP2VUgR9p5jZjGR9zexjggqu2sQ8/eZ8Guce7BzHrJnJQK1BXUFiiG58BZyjtuINlKEJp26T7DfUsBav1MNB+Q6aKR7FLfB1IzLaCAtNjLlSSqKLgkedXQJQbmBV4CxqyRXKhTn9ybLmPybRhOBA47gfRtavaLfed/BMqSRQjCcwCTwN36H6gGkjIsbf7mG1FUOuA8BVeGqEKoAssN9RYggYKTwEfwrqrUftrVMXZVOxeL4Oy47cOZniA/L3pR8P2MLSkxhI03bgBbY1G9Q1bLCIngumWY3Zb6WFQbuJ5271IH0/zlCVa3PZNqC6NUjhQ9IDScTnfQtZxZpolvEnsbmqJFFOqP/cwyHr5TxRrXW1135HWxuK2nvjhcXNy9eXJ1TOlXygdDYo1hVwF9mjtpYUDo0vke5d8Ukx/p9nDQJt47JVr2sT85esUBORSdu7km4FQfcoP6L+SmIwZutcfE9aS7zfyK154My6TmIwxwmmPw7rWc0kH9m/ykCOByZhCSxZIvs1bUiJbNy6ByZh6MUWif8ckSL5p9mNSZtFJkRice7IxuFwpmaZT/AN0/z3AXPJU5saeJ787EG8+bVkMJHuQZGhlcLlSMhX4n9mZWor2oNw7aGNmvvaat5gVcOgK9zhx9+BB5Z+b41bKO3aFWcNogPu4Wq7hNnMv6i5Oc0sInOKWcIPvbk4VmgvuH0AewBOVG2FfFxgeFpwTAdwHXtQYFytQFpEAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\int_{0}^{2 \\pi} \\sin{\\left (m x \\right )} \\sin{\\left (n x \\right )}\\, dx$$"
      ],
      "text/plain": [
       "2⋅π                     \n",
       " ⌠                      \n",
       " ⎮  sin(m⋅x)⋅sin(n⋅x) dx\n",
       " ⌡                      \n",
       " 0                      "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# m and n are symbols defined as integers\n",
    "sinIntegral = sp.Integral(sp.sin(n*x)*sp.sin(m*x),(x,0,2*sp.pi))\n",
    "sinIntegral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOAAAABwCAMAAAAwlYyQAAAAOVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACXHtMAAAAEnRSTlMAmc0y3e8iEESru4l2VKWLZpM8SRfxAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHdklEQVR4Ae1ciZKDKBBVURHFHZf//9jtbuQQEoXI7KgVqybx6OvRgBkfdlXZrW7s7hN3WNuxJ+IymJjqerP/xG/eqmHFxetpmfKSOQg5lmmVD5y/cRxYEspG2MxVxducdLKOD7Z93vhLPP3CeU4knpetJa6UydnSodSUM+HUtnE8Bx/tvnLe1J+YCixNimChpVHg56w4fqVtglTSZPelXjhnY5sRijUfWGqUjVHvDQo6aupWDuAL52KQMjUQTy6w1KlpvcgV9YhBJVtlousQIp+knGCmGppOzlOzDp1hEn0vpZDVMMnRTGReJP7uC+ccxkrryTh73sl4N7SkLJ5eQ2UWcawcnRE4BhnNSyO2S1fPXC1abKpGPMeVBHCzHQgwyL3NjLIXzifoSfhnttrYs3NPkiUfIHnLB9hSH2A4WTWuydlctThf6dAXZeJ8/d3r7uM7R1O9m/GcPTMppllyAMPcvtbfnsUMsnXQYuYbO6ArDrnDfqknLUr1VndzFDtfaKQIi8bYO2qq0JIDWOnROWdNMgBwXm+EeD/dTOwzZa2mpHZHAztyrmdQhmNg3bQ94ZJqLmy/A0sewJGMLVm3CVCxGYQevgGof8JTBx5UDynV29R6m8146HxeR6eZskBX24MhnmXJA6jvkFn/WVDP6ygW+kWzAdjh2NT9FIz23nyxQvW/QucGmAEKsmQPe7yZ9319tx9Y8gBWLcTAs/6zoHZnpNKg23Ztd3TXU9fVQwYu7kcFqhvnrmvam722Bz+0FjsuHSp/b2vJB4j3MzeqfZ3X+0wohXeCvp4mgffBEY5togYafJI64VwfRUU3U+cczenN3uy1vV4s1oMRCb63MHyAgeBfHlp8MMbPxaHUP+cMXF37ohks12xfgOXa8m8sfTP4N+1ezqtSP+WMXdHSt4teMSs5MX0zmNNaV5R9fAZ/fuwznCu2//mYHp/BL8D9TlKOXdr38/nVcxksyC59DmFf8xzAcuzSfpQnrp4DWI58OQFhX/Xcj+1bAHz/zFmqdfMeBvrNVZBd8s3SfiKXFOlFJ3a66DSxWfaL7N//FCjGLkVxxVxSJJJ24j1AoHOqqa/k3mM7AliCXQqDTeWSQr342APIG0caNBoWUA/CMAGxLlz0uQl4dv0xu1QFzlO5pFdBbc95ALcX8IjBY2nH+MXXCWAZdikynsglRXrhiV2ANZC1u2TTJoMwFW3Ilyx2KYwr5pIiibQTewA5gttdBkNjsAi7FEWbyCVFeuGJPYASuyc5CrXMcUF2yZhcv5O5pEAvOlRmzUB0BRbO4HIC8Z73KsouBf6TuaRALzzs1W4XDMXvdzwoIJefvE3qiNO/OfrWrgO6OZA34c9q9z7+RutGp1u7lvJGQWeEWj+8g04798CMZrqs6Pzw/gnLs5u9f4Uum5j0wHj77BcKcLXVw+8SldTrUNOTfjdJeKnAWzR1t+hT4h318tEU0XvKLOpgMfU9Ybmo2ck+enl2CQbh0fJL1xrx3pdditvk/z6z99DpOJabky+PA3jEJoWAS7JLxdikMEiviyawSaE2Pfgt8e5SVYxNCkN0AFPYpFCbAJZgl8qxSWGIDiBeOWKTQu0NN3GGXSrHJoUhbgAeskmhNgIsxC6VYpPCEDcAD9mkUHuTwXPsUvRmUujrw2Ol/rWax2ySFV13aAyWYZdKsUlhiH4Gj9mkULscu1SMTQpD9AEesUmhbkl2qRSbFMZY+QCji0848QV49yx+V/zePYPfMfjN4MVb4NtFL56gw/DOrfg9NP/3At8u+vc5OBfBN4Pn2u/vtf0MBhXXPgiOjVSzI02THRaZSbOzL+XPotuKa/t60VVNpPLOVoiJJMITfn2j8Fq5Yy+DQcW1PB98fbdCVzXL0/1VaQ9gUHEtz+1yA4Cazc4qHGcbYTarGTIy2LO9VzKs6ZM7LoNhOasUw6Zk3CJUKwRSqaPQxeJQu68XKB5n68ktrWpn+GiGSqhu6fX6ACaXRVJ9Mi2d4jZHxgF8UbvtyJBHu5hiRyM+mutpGfGA1YqGDj7WenJQQ84UeRpxnTG+RsLh2Ti9nOGkj7zmXfcB0jDKmds82sUCpBJ4tMpPr9fHe4GpJ0clw2iepZXu2J2pASqQcdJ5AI6kHcD8LmoLjkG4FiAtkcZUMVUvsAnIqHnjZ4bMzhLWdQx0U6Hx2kJZKEinJ30Uct51BzC/cJxPu1iA1OUQ4OKWN9gXYiBNddXNa4kuAsinRsEM5UnnATiS9gCGtduOVDcl4xAgoIICi6imM4jHtFmAdctlNY1cv4uBAHtc8LiogaHKb2ydmd/BDQ2fnMJxPu0CweIYdgB17cYKELlXmpiCWRWw6FcZEOBMuKAfU6VHki4Lc3QlfoPabQl+vJJxWMIwAMhoqqkhRab/wh4OUY1FN8aMs2wFk6qTTnCcIeKv5dpWXEsxYkrGgawQ8DYl0jGCzaOiqZKLWk58U09O4ttC9KFl+TAtcsHqq9UqneI2S2bv5aUsQ1cVlo9fFCt+a/66TEqbh79dAO9a59QtvkxecgIZ9I03R+UGsv8BiSpQOmIIf+8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\begin{cases} 0 & \\text{for}\\: m = 0 \\wedge n = 0 \\\\- \\pi & \\text{for}\\: m = - n \\\\\\pi & \\text{for}\\: m = n \\\\0 & \\text{otherwise} \\end{cases}$$"
      ],
      "text/plain": [
       "⎧0   for m = 0 ∧ n = 0\n",
       "⎪                     \n",
       "⎪-π     for m = -n    \n",
       "⎨                     \n",
       "⎪π       for m = n    \n",
       "⎪                     \n",
       "⎩0       otherwise    "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sinIntegral.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Because the integral is zero for $m \\neq n$, then the series on the LHS:\n",
    "\n",
    "$$\n",
    "\\sum_{k=0}^N \\sin{\\left (\\frac{\\pi x}{L} k \\right )} \\frac{d a_k{\\left (t \\right )}}{d t}\n",
    "$$\n",
    "\n",
    "is a linear system.  This is also true for the RHS:\n",
    "\n",
    "$$\n",
    "- \\sum_{k=0}^N \\frac{\\pi^{2} k^{2}}{L^{2}} a_k{\\left (t \\right )} \\sin{\\left (\\frac{\\pi x}{L} k \\right )}\n",
    "$$\n",
    "\n",
    "The principle of superposition permits us to split this summation into N independent ordinary differential equations, solve each one, and then sum the solutions to produce the answer to the original PDE.  To continue developing a solution for these N independent ODEs it is necessary to analyze the amplitude ODE and define a differencing scheme.  This will be illustrated using SymPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMcAAAAxBAMAAACc4KDxAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiKZiXbvRGa7VKtFbb1tAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD90lEQVRYCbVYT0gUYRR/s/9m/7S7A9Uxd9lI6JBKEpWHWIKgDtGeK2pLKA8hm6lXN6IOIWSnDh3cIKgOope6RGkdwvCyRAcPiYt08CChlBZmTu/7MzPfzHwzO5p9h2/e+73f+733dnZnPgWwr3uGG8kb1s5fT3HJSG/HzouLihnipOVFzEnFBGYPftAM0IdlUHYRw6OIManBta5qKZM3PG8WYTxv4OZXhJDkK1yKrMojTvTmCiLbLFKLrjvlpH50H7kTPkXopNJUBNOkQ7J8WahfhYXCwUKh6nFP6KRMyb2HqxzzZQHk9iDPPsmsoMYmFQCbOcg9fxbAJxWJrMgITxGLkEk9V6pshPxYyGkZwI0UiXZvVvGCy1aETspw1/4MnnLMmxW93vOQkegkpoZRhMbppGbIZkQOzb3jgDfrRSm0zEjyIjROJ7VpE+eBrutrCV3HfNqKnIVE5THEpvCKK8oufOeTCHFbGB21f/bC3RJDrVZFlrK/E9fhYmwFzK+gFe8rFI4XCm0IyOOUeRYG4kWW49MKISSmIMeZjG/sfBIWjz8yYNt1KqMx392KeuOKRc1VYB5AIsKLsDh0WRmWFWmEueNuZQAWLGKuDu/TJYmIUYTGFePJYSWiFZvI8kncrRyBc/x2ITFUTZ1QQSLCi7B4ssOmzp2clq0xk7UqqnTBZR5DhrLYf+k+SER4ERZXy6cbTE7c70CiwnxZK/N8SjNDIsKLcJE3TyQf2AxkeBFZK8dMdW6E3CJ7RU72q+h52UxF+Ug/W7XspDUTOTNcd6ZIfK4SovKvXIRmIsPJkaIryQVwlSxpKFlJOjOaiUxHOvin71IWAK4yriF2cu4LuYirmcgBWCyKfLnNVW6R6KiuO0nBRJxZTp+rtMJsoO+JM9301QJdZRNwG8qf+PnfbnhnkfTKbWWCSeL7ZqvL2Ys8H2JrTuJ/8MOL0/+oqpJXXWdnw0cmVxMejT68bYSss/0k5OoBflTbqCGc7ZcgW/Qvwr8WWy4jnO1nIN6j+QlEfvlF5bG3elsdwkHP9qgRG6U97C7J9eTodwqbZ3s5SUDHJivEc73UBIrLTP6kkORg5aIyoJygL8CrHmEprK5T2DjbSzkimKynaFtHRbCZzU651tm+GT+Gj2gNkksbi82YQjybJ451thdCUnMMgNyU9LI06gGOk9sonO09aCbcwGMr3pRY3kQCGEM1JLGzfQA2pOsAKXyGhklvgVc7MpXAbKD/ksCbkqsFzwHYQPLL4AlzhDpZhvGtdJYkD4kyyQy2GoSWWIU++BwsgbDUH/hHwURgfrr1Gq7uNehVqoGTILwMMKQF5of5+1S7+DpwDrS0b35r35Tz/wIMDUSHaabMSQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{- a^{i+1}_k + a^{i}_k}{\\Delta t} = \\frac{\\pi^{2} a^{i}_k}{L^{2}} k^{2}$$"
      ],
      "text/plain": [
       "                          2           2\n",
       "-a_k__{i+1} + a_k__{i}   π ⋅a_k__{i}⋅k \n",
       "────────────────────── = ──────────────\n",
       "       \\Delta t                 2      \n",
       "                               L       "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ai, aip1 = sp.symbols('a^{i}_k, a^{i+1}_k')\n",
    "dt = sp.Symbol(r'\\Delta t')\n",
    "\n",
    "differenceEquation = sp.Eq((ai-aip1)/dt,((sp.pi**2*k**2*ai)/L**2))\n",
    "differenceEquation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKcAAAAyBAMAAADLk9YvAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD9ElEQVRYCa1XTWgTQRR+u9mkiel2g4cqKCYiKj01Kr2IaBREQbAq4kEUo4gggs1FURCbi1BEbQ7iwR+siCetxpOipRY8qCcLXrw1oggWpK0V/Kka572Z2Z0wu9k1+A7zfr9vZ19m32YBmMzP4/o/5GpZslh7StL01bZv1Ddo782I+GDIRvtrvgS+wSsHeTj10zftBXtm0DaOepFgKzXHc2nCBNdZF79gMhlSJgg28FY6Y8IPUDE4gI2KhZQJtOiVUw1gE+GtMJFjZgKXcMnyHygbQpqHOLvzbQ/u1sIpAQRbCKlVBhubOhiFMippmpFhU9e1TPpeRw6w0ARr1CEvZa5dJR0PsOw2u7LvTk/LYk9jx1lTja9eqAtGpOMCzEp6DIOReposs0p7FqxpoJ0g8Cj0V1Crkqi0/UA/Eim2lDW1Yo4lS2Sy5TDcLEhb6kQhRfeik27He22UC2tQ+vKx6hIlMSh2rQL4I6eTLtUfxb46yUxqeFwhPS5sFZCoYVAjTS0sYjxczDyv4QDjThHdZRTTSGFejRKhyxtZwQHteA2bX0gnzS6Q1U21lbNKvIADnDLzNsEWjIlJcjOHDskFU1pN9ZmB8xkV0Mu8tpcDtzDmDFFmokAKl0VdrtnMOFCvizQHrGBevF6fxmA7rXCEClKr1+8ko+kyxY7CLK9QAM/h0iMJS/5Gy6QHATZXxDVk1k+bKy/d+FThGQ9gfO+4TFSUeIb55QW0jf0QG0KjqZyCro4Sr1AAyZl3xriLS+8GWPyQ3NgMsMPbsYucxUdQXqDt8LPPV+YPpTNUAgogJhrCE3CuBp18isWHIMv2cFgkglRbNSFSCiAxel+tP5kD6KH7yeZwtBv6Y6qWs7ffuCN2qgCyBXW+mOz2IfYNcdky3EtWrCLafiLmcjbjFHhaAUwwdM4FddMm6Ui11+xjJpj5s1U322CIufwR4gKuAJ6AUxJRAOsP4aYKTBmjK699gPa3GwMaIObyMKQFXAEMQ8d60RSAedNEqjymjnuGKUOL/ezXY/+5zNI6QPw1Ea9/ZDgxWUbVKE6RfDmXG5I6QJ9Sk9YOanMDsDdHrpzLDTkdoJPebytyBhXZXUBPzmU1A6ADdNJXMKrv9DXxuHO5gVUH6KQNAOn8QsObyzIcoKORWvTv1JvLAWQyHI3UpNe5N5clOkBHI6Ux657tACovHI0UXzp2xUOFWNFIe2ts6IQwKelopDh1PiuoEDMaaR87+AFTxo9fkDb9kDBG6k9H6lU/uH9MTJL4mH+6tWh/jnBp/o5ujUNDTbGfgIn88tPyrQQM+U06VWsF7o/p3Cni9lzGv+Lfo8l9LtWif/h1m1/o+jjm/wKGGBfIfgfECAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left\\{a^{i}_k - \\frac{\\Delta t}{L^{2}} \\pi^{2} a^{i}_k k^{2}\\right\\}$$"
      ],
      "text/plain": [
       "⎧            2                    2⎫\n",
       "⎪           π ⋅\\Delta t⋅a_k__{i}⋅k ⎪\n",
       "⎨a_k__{i} - ───────────────────────⎬\n",
       "⎪                       2          ⎪\n",
       "⎩                      L           ⎭"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "odeSolution = sp.solveset(differenceEquation,aip1)\n",
    "odeSolution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The above solution results in the following difference scheme:\n",
    "\n",
    "$$\n",
    "a^{i+1}_k = a^{i}_k \\left( 1 - \\frac{\\pi^{2} k^{2}}{L^{2}} \\Delta t \\right)\n",
    "$$\n",
    "\n",
    "The timestep, $dt$, should be chosen small enough such that the $a_k$ decay at each timestep."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Implementation of the Spectral Method\n",
    "----\n",
    "\n",
    "An annotated implementation of the spectral method is developed in the next section.  Unlike the development above, we use the full Fourier series for the basis functions.  This is a basic implementation that could be improved with the addition of helper and visualization functions.  `SciPy` provides $\\sin$ and $\\cos$ transforms for other boundary conditions where Fourier may be inappropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The following quantities are defined:\n",
    "\n",
    "* `numPoints`: the number of points in the grid\n",
    "* `L`: the total length of the domain\n",
    "* `dt`:  the numerical timestep\n",
    "* `steps`: the number of timesteps to compute\n",
    "* `x`: a vector containing the x-position for each grid point in the domain from $0$ to $L$ with `numPoints` entries in the vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "numPoints = 20\n",
    "L = 1.0\n",
    "dt = 0.0001\n",
    "steps = 1000\n",
    "# we have access to np.pi for $\\pi$\n",
    "x = np.linspace(0,L,numPoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The following vectors are defined:\n",
    "\n",
    "* `c_new` will hold the $c(x,t)$ values at the start and end of the numerical computation\n",
    "* `a_old` will hold the amplitudes at the current timestep (i.e. the $a^{i}_k$)\n",
    "* `a_new` will hold the amplitudes at the next timestep (i.e. the $a^{i+1}_k$)\n",
    "\n",
    "They are initalized to zero at the start of the calculation.  Each of these vectors is the same shape as `x` to hold the results of the FFT and IFFT.  They are declared as `complex128` to accomodate different initial conditions where the FFT would produce complex valued results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "c_new = np.zeros((numPoints), dtype='complex128')\n",
    "a_old = np.zeros((numPoints), dtype='complex128')\n",
    "a_new = np.zeros((numPoints), dtype='complex128')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In previous example calculations a set of basis vectors was created to illustrate Fourier transforms graphically and the terms in a Fourier series.  That is not part of this calculation, however this could be easily added if visualizing the results in more detail is desireable.\n",
    "\n",
    "The difference equation requires values for the square of the Fourier numbers, $k^2$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "a^{i+1}_k = a^{i}_k \\left( 1 - \\frac{\\pi^{2} k^{2}}{L^{2}} \\Delta t \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To ensure appropriate matching between the $a_k(t)$ and $k$ the helper function `fftfreq()` is used.  Using `fftfreq` with `fft` and `ifft` ensures that the bookeeping of amplitudes and $k$ values is done correctly.  In the equation for the time rate of change in the amplitudes there is a constant $k^2$ term - computing this first avoids having to repeatedly compute the quantity each time step.\n",
    "\n",
    "`NumPy` arrays are called by reference so it is necessary to perform an element-by-element \"deep copy\" of the data from one array into another array when building the initial condition. The helper function `np.copyto` provides this capability.  `a_new` is then filled with the amplitudes corresponding to the initial condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "k = np.fft.fftfreq(numPoints, d=L/(numPoints-1))\n",
    "k2 = k**2\n",
    "initialCondition = np.sin(np.pi*x/L)\n",
    "\n",
    "# create an initial condition (this could be a simple function like x**2)\n",
    "np.copyto(c_new, initialCondition)\n",
    "# transform it (dft or sin transform)\n",
    "np.copyto(a_new,np.fft.fft(c_new))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Instabilities will occur if the amplitudes do not decay at each timestep.  The problem is that the condition depends on the wavenumber - so a suitable $\\Delta t$ must be chosen that satisfies the most restrictive condition for the largest wavenumber.  Using a Boolean expression it is possible to check to see if all of the wavenumbers result in a factor less than one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True,  True,  True,  True,  True,  True,  True,  True,\n",
       "        True,  True])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(dt*np.pi**2*k2)/L**2 < 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If any the results are `False` then the numerical calculation will not converge, if they are all `True` then it is possible to complete the numerical calculation.  The next code block performs the numerical iterations of the amplitudes.  First the pointers to `a_new` and `a_old` are swapped and then `a_new` is filled with the new values based on the `a_old` values.  This sequence of operations is performed for the number of timesteps given in `steps`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(steps):\n",
    "    # swap pointers\n",
    "    a_new, a_old = a_old, a_new\n",
    "    # find new amplitudes\n",
    "    np.copyto(a_new, a_old*(1-(dt*np.pi**2*k2)/L**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When the requested number of steps have been computed, we use the inverse Fourier transform to compute the concentration profile and store those results in `c_new`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# inverse transform it\n",
    "np.copyto(c_new, np.fft.ifft(a_new))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "After the computation, the concentration profile are displayed with a helper function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def makePlot():\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8]) # left, bottom, width, height (range 0 to 1)\n",
    "    axes.plot(x, c_new.real, 'r')\n",
    "    # Setting the y-limit cleans up the plot.\n",
    "    axes.set_ylim([0.0,1.0])\n",
    "    axes.set_xlabel('Distance $x$')\n",
    "    axes.set_ylabel('Concentration $c(x,t)$')\n",
    "    axes.set_title('Concentration Profile solved by Spectral Method')\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAElCAYAAAAhjw8JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3XmYFNXVx/HvAVwRAQVFAQXFDXGDUcEtqGjQl6hxi+KuwWhckohxj7txi6JxiUvcdzEuaDQYDYiCIIMLioiiouAKLoggynLeP06NtO0M0z013T098/s8zzzTXVVdfaq7uk7VvXXvNXdHREQkjWalDkBERMqfkomIiKSmZCIiIqkpmYiISGpKJiIikpqSiYiIpKZkInViZmuZ2bdm1rzUsSyNma1gZo+b2WwzG2pmB5nZ0xnz3cy6FTGekWb22wKs93Yzu7CGeYeb2Qv1/Z7lxMzONbO7S/C+08ysXz2tqyD7Tn1RMqkDMxtoZpXJwfQTM3vKzLYrdVzVMbMuyQGzRcr1/ORH4e4fuvtK7r4ofZQ/ey83s7nJ5/uRmV2ZImntC6wOrOru+7n7Pe6+az2G26SY2cZm9rSZfWVmX5vZBDPbvcDvWdCDqJn1Tfa5h7Omb5ZMH5njempM6E2BkkmezOwk4Crgr8RBai3gemDPUsaVRtpEUyCbuftKwM7AQGBQ9gI5xr028La7L6zn+Jqqx4H/Evv+asCJwDelDKie9t+ZwDZmtmrGtMOAt+th3U2Du+svxz+gNfAtsN9SllmOSDYfJ39XAcsl8/oCM4DBwOfAJ8ARGa9dAbgC+ACYDbwArJDM6w2MAb4GXgP6ZrxuJHABMBqYAzwNtEvmfQh4Eve3QB/g8GTZIcCXwIXAusD/gC+AWcA9QJtkHXcBi4HvknWcAnRJ1tsiWWZNYFiyvqnAoIz4zgUeBO5M4psEVCzlM3SgW8bzocC1yeNpwKnAROB7oAWwUfIZfJ2se49k2fOAH4AFSdxHJdv+QnXvlXx3f0s+s8+AG6o+/2pi7AY8l3xPs4AHMuZtA4xP5o0Htsn6rn6bvNfXQI+Mee2Tz3i15PkA4NVkuTHAphnLbgG8nHyeDwD3AxfWEGvV931NEtNbwM7JvP2ACVnLDwYerWY97ZLPq00N79OX2L/PSD6TacBBWb+NGj9f4oTsVSI5vQv0By4CFgHzk++waj9w4DjgHeD9ZNrVwPTk9ROA7bP2wbtrifsG4LhkWvNk2tnAyIxlNySS6ZfAFGD/ZPrRxH72QxLn4xn768nE/jo7+a6Wz1jfIOL38iXx+1kzY94uyXc1G7iW2N9+W+rjYI2/21IHUE5/yc69kOQAWsMy5wNjibO29slB4IJkXt/k9ecDywC7A/OAtsn864iDTcdkZ94m+QF2JA7yuxNXk7skz9snrxuZ/PjWJxLSSOCSZF4XMg76ybTDkzhOIA7GKxAHx12S92sPjAKuynjNNKBfxvOfrDfZ0a8Hlgc2J870qg5Y5xIHg92T7boYGLuUzzDzAN8d+BQ4KiOOV4HOSdzLJD/GM4BlgZ2IA+wGGe99d9a215RMrkp+0KsArYiz8ItriPE+4Mzk+1ge2C6ZvgrwFXBI8tkemDxfNeO7+m3y+Fbgoox1Hgf8J3nckzjh2Dr5zA5Ltn25ZDs/AP6UbP++xIFsaclkYcbyvyEOUKsk6/sS2Chj+VeAfapZjxEH7yeAvYDVs+b3Td7nymS9vwDmZnwXNX6+wFZJTLskn2lHYMPszyzre/tvsq6qE66DgVWTz30wsd8sX91+UE3cM4jf27hk2u7AcCLxj0ymtSSS1RHJe/QkkubGyfzbs7+D5Dt7iTjZWgWYDByTzNspeX3P5PO6BhiVzGtHJMV9k+/sT8lnq2TSGP6Ag4BPa1nmXWD3jOe/BKYlj/sSZ56ZB/bPiauOZsm8zapZ56nAXVnThgOHJY9HAmdlzPs9Sw5KXag+mXxYy3bsBbyS8XwaNSQT4sC+CGiVMf9i4Pbk8bnAMxnzugPfLeW9PfkhfZV8nhcCzTLiODJj2e2Jg0azjGn3AedmvHetyYQ4UM4F1s2Y14fkrLeaGO8EbgI6ZU0/BHgpa9qLwOEZ31VVMukHvJex3Gjg0OTxP0hOQjLmTyEO0DsQV72WMW8MS08m2cu/BByS8V4XJY83Tj735WpYVyfiLPld4mp1FLBexv69EGiZsfyDwF9q+3yBG4EhNbznj59Z1ve2Uy378Fckv6fs/SBrub7AjOTxO8AGxJXeQfw0mfwGeD7rtTcC5ySPb8/+Doj99eCM55cBNySPbwEuy5i3EnFS0AU4lIwTruTzm5H9OTSkP9WZ5OcLoF0tZbRrEmeNVT5Ipv24Dv9p+f08YidqR5zhvlvNOtcG9ksqPL82s6+B7YA1Mpb5tJp1Ls30zCdmtpqZ3Z9UeH8D3J3ElIs1gS/dfU7GtA+Is8ua4lu+ls+xp7u3dfd13f0sd19cQ+xrAtOz5me/dy7aAysCEzI+4/8k06tzCvEDf8nMJpnZkRnxfJC1bE3x/A9Ywcy2NrO1iSu6R5J5awODs77zzsn61wQ+8uQok/EeS1Pd8lX75R3AQDMzIhk+6O7fV7cSd5/h7se7+7pJjHOJxFrlK3efW8371Pb5dqb6fX9psvfhwWY2Oblz72uiWDrXfbjKXcDxwI4s+S6qrA1snfWdHAR0qGWdNf02f7KvuPu3xDGmYzJvesY8J2t7Gxolk/y8SBTX7LWUZT4mdroqayXTajMrWfe61cybTlyZtMn4a+nul+SwXs9x+sXJtE3dfWWiyMByWA/E9q1iZq0ypq0FfJRDfHWRGcvHQGczy9yX6/Les4grw40zPuPWHjcB/DwA90/dfZC7rwn8Drg+ucU4+/uvMZ4kAT5IFIUNBJ7ISMjTiauFzO98RXe/j6hr65gc/DPfY2mqW/7jJI6xRFn/9kkcd9Wyrqr4pxNFsz0yJrc1s5bVvE9tn+90qt/3IYd92My2J67g9yeKjdsQxWZWw2trchdxZf+ku8/LmjcdeC7rO1nJ3Y+tJc6a/GRfST63VYl95RMiwVbNs8znDZGSSR7cfTZRIXedme1lZiua2TJmtpuZXZYsdh9wlpm1N7N2yfK13t+eHFhuBa40szXNrLmZ9TGz5ZLX/8rMfplMXz65nbFTDmHPJIoj1qlluVZExeHXZtYR+HPW/M9qWkdyUBkDXJzEtilR2X1PDvGlNY44Oz4l+S76Ar8iiilylnz+NwNDzGw1ADPraGa/rG55M9sv4/P/ijiQLAKeBNZPbh9vYWa/IYr1nqjhre8lik8OSh5XuRk4JrlqMTNraWb/lyTsF4nipBOT99ibqHNYmtWS5Zcxs/2ImxaezJh/J1F8tdDdq22TYmZtzew8M+tmZs2S/ftIoo4w03lmtmxygB8ADM3h870FOMLMdk7W3dHMNkzm1bjvZWiVfCYzgRZmdjawci2v+Rl3f58oSjyzmtlPEN/tIcnnuIyZbWlmG+URZ6Z7iW3ePPmd/5Wos5kG/BvY2Mz2Tq7gT6T2K6CSUjLJk7tfCZwEnEXsuNOJy+JHk0UuBCqJuzdeJ+64yfXe85OT14wnKkUvJeoCphN3upyR8Z5/JofvLzm7uggYnVya965h0fOIisDZxI78cNb8i4kk+bWZnVzN6w8kyno/JooHznH3/9YWX1ru/gOwB7AbcfZ7PVHv8FYdVncqUZk/Ninqe4YoP6/OlsA4M/uWqFT+g7u/7+5fEAfQwUSRxSnAAHefVUP8VclwTeCpjOmVxJ0+1xLJaipR91G1zXsnz78iklH295VtHLAe8RldBOybxFrlLuIKY2lXJT8Q3/EzRJ3WG8QddYdnLPNpEtPHxMnEMRnfRY2fr7u/RFRsDyH2wedYctZ+NbCvRduWv9cQ23Di83ubKDqaTx2Lhdz9BXf/WWlCctW4K3BAsn2fEr/R5ZJFbgG6J7+RR7NfX836niXqk/5FXImsm6ybZH/ZD7iE2I/WI+rUGiz7aTGqiDRFZrYCcTNIT3d/p47r6EtUcudyxSyNjK5MRATgWGB8XROJSFGSiZndamafm9kbNcw3M/u7mU01s4lm1rMYcYlIdJUD/IEomhOpk6IUc5nZDkTl7p3u3qOa+bsTDeh2JxppXe3uWxc8MBERqRdFuTJx91FEhXJN9iQSjSe3KbYxszWWsryIiDQgDaWDv4789M6LGcm0T7IXNLOjiX5waNmyZa8NN9wwexEREaknEyZMmOXuNTXe/VFDSSbVNSyqtvzN3W8iurGgoqLCKysrCxmXiEiTZma19a4ANJy7uWbw09adncit1biIiDQADSWZDAMOTe7q6g3MdvefFXGJiEjDVJRiLjO7j+iZs52ZzQDOIbpVxt1vILp12J1oHTuPaAkrIiJloijJxN0PrGW+E2M5iIhIGWooxVwiIlLGlExERCQ1JRMREUlNyURERFJTMhERkdSUTEREJDUlExERSU3JREREUlMyERGR1JRMREQkNSUTERFJTclERERSUzIREZHUlExERCQ1JRMREUlNyURERFJTMhERkdSUTEREJDUlExERSU3JREREUlMyERGR1JRMREQkNSUTERFJTclERERSUzIREZHUlExERCQ1JRMREUlNyURERFJTMhERkdSUTEREJDUlExERSU3JREREUlMyERGR1JRMREQkNSUTERFJTclERERSUzIREZHUipZMzKy/mU0xs6lmdlo189cysxFm9oqZTTSz3YsVm4iIpFOUZGJmzYHrgN2A7sCBZtY9a7GzgAfdfQvgAOD6YsQmIiLpFevKZCtgqru/5+4/APcDe2Yt48DKyePWwMdFik1ERFIqVjLpCEzPeD4jmZbpXOBgM5sBPAmcUN2KzOxoM6s0s8qZM2cWIlYREclTsZKJVTPNs54fCNzu7p2A3YG7zOxn8bn7Te5e4e4V7du3L0CoIiKSrxZFep8ZQOeM5534eTHWUUB/AHd/0cyWB9oBnxclQpG6WrwY3nwTXnoJ5s9Pt65VVoFtt4XOnWtfVqQBKVYyGQ+sZ2ZdgY+ICvaBWct8COwM3G5mGwHLAyrHkoZn4UJ49VUYNSr+nn8evvyyft+ja1fYYQfYfvv4360bWHUX+CINQ1GSibsvNLPjgeFAc+BWd59kZucDle4+DBgM3GxmfyKKwA539+yiMJHimz8fxo9fkjzGjIFvv4153brBXnvFQX/bbaF163TvNWNGJKdRo+Df/4Y77ojpHTpEUqn623hjaKZmYtJwWDkfrysqKryysrLUYUhj8+238OKLS5LHuHHw/fcxr0ePJQf07beHNdcsXBzu8NZbS+J47jn46KOY17btkquWHXaALbaAFsUqaJCmxMwmuHtFrcspmYgQyePhh+OgPWECLFoUZ/49ey45YG+3Hay6aulidIdp05YUrY0aBe+8E/Natowrox12gIEDo5hMpB4omYjUZvHiKEq69FIYPRqWXRa23npJ8ujTB1q1KnWUS/fJJ0sSy6hR8Prr0Lw57L8/nHIKbL55qSOUMqdkIlKTH36Ae++Fyy+Pu7DWXhsGD4Yjj4wz/HI2YwZcfTXceCPMmQO77gqnngo77qgKfKmTXJOJavCk6ZgzB668EtZdF444IuoY7r47iopOOKH8EwlAp06RJD/8EP76V3jtNdh5Z9hqK3jooSi+EykAJRNp/D77DM48E9ZaK65AunWDp56K23sPOgiWWabUEda/Nm3g9NOjjuXGG+Hrr2G//WDDDeN52vYwIlmUTKTxmjoVjjkmirEuvjjO0MeNgxEjoH//plHss/zycPTRcVfY0KGRZI45Brp0ic/k669LHaE0Ekom0vhMmBAV0BtsALfdBoceGgfThx6K4p6mqHlz2HffaKX/v/9FxfwZZ0RL+5NPXnLLsUgdKZlI4+AOTz8dVx8VFTB8eNzNNG0a3HQTrL9+qSNsGMyiMv4//4FXXoE99oCrropbiY84Im5IEKkDJRMpf6NHQ69e8MtfwuTJcNllMH16FOOssUapo2u4Nt8c7rknbkD43e/ggQeiZf0++8TnJ5IHJRMpX3PmwPHHR0vwL76AW26B99+HP/8ZVl659tdL6NoVrrkGPvgAzj47rlo23hiuvz7a4ojkQMlEytOTTy454J1wAkyaFO1Elluu1JGVr/bt4bzz4I03ovHmccfBL34BU6aUOjIpA0omUl5mzYKDD4b/+z9YaaUo4rr66ngs9aNr16h/uu22SCybbRZtVhYsKHVk0oApmUh5cIf77oONNoIHH4zimFdeiS5PpP6ZweGHRx3Ur34V7XS23DLulBOphpKJNHzTp8cBbeBAWGcdePnlKI5RkVbhdegQ7VMeeQQ+/zxurT7lFJg3r9SRSQOTdzIxs5Zm1rwQwYj8xOLF8I9/RN3IiBHRFcqYMdENvBTXXnvFbcNHHhndtWy6aXwnIolak4mZNTOzgWb2bzP7HHgL+MTMJpnZ5Wa2XuHDlCZnyhTo2xd+//uoDH7jDfjTn6LxnZRGmzZw883R6BFgp52idb1a0Qu5XZmMANYFTgc6uHtnd18N2B4YC1xiZgcXMEZpShYsiMrezTaL7tRvvTUqgzU+R8Ox444wcWK0nL/lFujeHR59tNRRSYnlkkz6ufsFwGx3//Gmc3f/0t3/5e77AA8ULEJpOiZMiEreM8+MOpLJk6NVdlPoQ6vcrLhiFHeNGxe3FP/619GFzWeflToyKZFak4m7V90P+Ej2PDPrnbWMSP7mzYtK3a23joPRww9HpW+HDqWOTGpTUQGVlXDhhfDYY3G33e23x9130qTkUmeyv5ldArQys42yKt9vKlxo0iRMnhwHpMsvj6uQyZPjLFfKxzLLxNXka69FkdcRR8Dee8Ps2aWOTIool2Ku0cCbQFvgSuAdM3vZzJ4AvitkcNLI3X9/FGvNmhX1IjffHJW8Up423DCGDr7iCnjiiegv7bXXSh2VFEkuxVwfufudwJ7uvpu7rwP0A84Bdip0gNII/fBDdIFy4IHR2eArr8Auu5Q6KqkPzZrBSSfByJHw3XfQu3cUe0mjl3M7E3cfnfH4S3ef4O5zCxOWNFrTp0d/T9deG7f6jhgBHTuWOiqpb9tuGycJ22wTxV6DBml0x0ZOLeCleP77X+jZMzplHDo0GiE2xiFzJay2WhRfnnEG/POfkVjee6/UUUmB1DmZmNkaZqb+LKR2ixfD+efHeCOrrw7jx8eof9L4NW8OF10Ejz8ewwP06hWPpdFJc2VyF/CWmf2tvoKRRuiLL2DAADjnHDjooGiXsMEGpY5Kim3AgOhTbZ11YnTHM86AhQtLHZXUozonE3fvB6wD3FZ/4UijMn58FGs9+2z0sXXnndCyZamjklLp2jWGDBg0KEbB3HVXNXJsRHJOJmZ2ldlPmyJ7mFT/YUlZc4/ksd120Xp99Gg45hi1ZBdYfnm46aa4w+vFF2GLLeCFF0odldSDfK5MvgWGmVlLADPb1cxG1/IaaWrmzoVDD40OGnfeObpIqagodVTS0Bx2WBR5tmwZHXoOGaJW82Uun1uDzwLuA0aa2QvAYOC0QgUmZWjKlOgS5Z574IILouHaqquWOippqDbdNLpi2WOPaJuy//7wzTeljkrqKJ9irp2BQcBcoD1wors/X6jApMw89FC0Zv/sMxg+HM46KxqwiSxN69bwr3/B3/4WA3BtuWX0Fi1lJ59f+5nAX9y9L7Av8ICZqQV8U7dgQZxV7rdfDGL18stqzS75MYPBg2OclG++iavbu+4qdVSSp3yKuXZy9xeSx68DuwEXFiowKQNffBFtR4YMgRNPhOeeg86dSx2VlKsddohW81ttFfVugwfDokWljkpylEuvwdXeguPunwA7L20ZacQmTYof/Zgxccvv1VfDssuWOiopdx06wDPPxMnJlVdG+xSN5FgWchpp0cxOMLO1Miea2bJAHzO7AzisINFJw/TEE9CnT4xD8txzcMghpY5IGpMWLeLk5Kaboo1S797w9tuljkpqkUsy6Q8sAu4zs4/N7E0zew94BzgQGOLutxcwRmko3OHSS+Pum/XXj0aJW29d6qiksRo0KJLJF1/Efvb006WOSJYily7o57v79e6+LbA2UbTV093XdvdB7v5qwaOU0vvuu7gCOe00+M1v4PnnoVOnUkcljd3228dJy1prwW67xRWL2qM0SPncGlwBLOvun7h73oWYZtbfzKaY2VQzq7Z9SjKq45tmNsnM7s33PaRAPv44uo2/557otO/ee2GFFUodlTQVXbpELwp77gl//GNcsXz/famjkiz53Bp8J1HcBYCZtTOzAbm8MBnq9zriDrDuwIFm1j1rmfWA04Ft3X1j4I95xCaFMn583Ps/eTI8+mh00Kf7LaTYVlop2jL95S9wyy3Ru8Lnn5c6KsmQTzKZ7+4/jm7j7rOA83N87VbAVHd/z91/AO4H9sxaZhBwnbt/laxfe0qp3XtvFDMsu2zctbVn9lcmUkTNmsVQBg88EO2ZttwSXlUpe0ORTzJ5z8x2y5qW672gHYHpGc9nJNMyrQ+sb2ajzWysmfXPIzapT4sXw+mnR5fxvXvH1ckmm5Q6KpGw//7ROeTixTGi47/+VeqIhPySyQnAxWZ2r5mdaGY3Au/m+NrqykWya9FaAOsBfYm7xP5pZm1+tiKzo82s0swqZ86cmXv0kptvvoG99oJLLoHf/S7uoGnXrtRRifxUz55xkrPZZjHQ2nnnRXKRksmnBfwnQC/gX0TfXBOBgTm+fAaQ2TS6E/BxNcs85u4L3P19YAqRXLLjuMndK9y9on379rmGL7l4770YWvXJJ+G66+CGG9QQURquDh1gxIjogfjcc+Muw7lzSx1Vk9Uin4XdfRGRTPK9rhwPrGdmXYGPgAP4eSJ6lLgiud3M2hHFXhowulhGjFgylO7TT8NO6nZNysByy8Ftt0UPxH/+M0ydCo89FrcSS1EVpVtXd18IHA8MByYDD7r7JDM738z2SBYbDnxhZm8CI4A/u/sXxYivyfvHP2LUuw4d4KWXlEikvJhFZ6NPPBFX11tuGbcSS1GZl3EDoIqKCq+srCx1GOXrhx/gD3+I4qwBA6Idycorlzoqkbp7663ooWHatNivjzyy1BGVPTOb4O61jnCXczGXmS0H7AN0yXydu+d6e7A0JJ98Et3Gjx4Np54ajRGbNy91VCLpbLhhjOB4wAFw1FEx+NZVV6nurwjyKeZ6jGgbspAYIKvqT8rNiy9Cr17R3fcDD8SdW0ok0li0bRs3kZx6ahTh7rhjnDxJQeVTAd/J3dX2o5y5R0+sJ5wQFZTDh6v9iDROzZvHSVLPnlHU1atXtKDfZptSR9Zo5XNlMsbMdOQpV/PnR59GxxwD/fqpIaI0DfvvD2PHwoorQt++UY9SxvXEDVk+yWQ7YELSWeNEM3vdzCYWKjCpRzNmREeNt9wSY7M//ngUBYg0BT16xMlTv35w7LFxUjV/fu2vk7zkU8yV3ZWKlINRo6Kifd48ePhh+PWvSx2RSPG1bRsnUeeeCxdeCK+/Ht2waBiFepNPC/gPgDbAr5K/Nsk0aYjc4ZpronfVNm2i/YgSiTRlzZvDBRfESdWbb0Y9yqhRpY6q0chnPJM/APcAqyV/d5vZCYUKTFL47js4/PAYR3u33SKRbLRRqaMSaRh+/ev4TbRtGydbf/+76lHqQT51JkcBW7v72e5+NtCb6DZeGpIPPoDttoM774zO7x59FFq3LnVUIg3LRhtFe5Tdd4+Gu4cdFidhUmf5JBMjY3Cs5LFGSWpI/vc/qKiI/okefxzOPjvGgBCRn2vdGh55JE667rorTsI+UMl9XeVzpLkNGGdm55rZucBY4JaCRCX5cYcrr4RddoH27ePOlQE5DYIp0rQ1axYnXY8/HidhvXrFSZnkLZ8K+CuBI4Evga+AI9z9qkIFJjmaOxcGDoTBg2McknHjYP31Sx2VSHkZMCBOwlZfPU7KrrhC9Sh5yrcL+gnAhALFIvl6++247ff11+Gvf4XTTtP47CJ1tf760cDxiCPg5JMjudx4o+occ1TrlYmZvZD8n2Nm32T8zTGzbwofovzM4sVw7bWw+eYwfXr0Q3T66UokImm1agVDh8bJ2dChMU7KiBGljqos1JpM3H275H8rd18546+Vu6u/8mKbPh1++cvoX+sXv4A33oD+6jJNpN6YxcnZ6NEx+NZOO8Gf/qS7vWqRTzuTS3OZJgXiHnecbLJJ9Pp7ww1xRbLmmqWOTKRx6t07etY+/vjoxr5q3HmpVj53c+1SzTR1sVIMM2fGkLqHHhr9DL32GvzudyrWEim0li2jJ4mnn4Y5c6BPn+iSZcGCUkfW4ORSZ3Ksmb0ObJB08Fj19z7weuFDbOKGDYsE8sQTcOml8NxzsO66pY5KpGnZZZcoUh44MNql9OkTXbLIj3K5MrmX6ItrGEv65foV0MvdDypgbE3bN9/EOAx77glrrBEjxp1yigaxEimVNm2iZ4mHHophgXv2hCFD4oYYyakCfra7T3P3A4FvgNWBtYEeZrZDoQNskkaOjLtI7rgDzjgj+hHS2CMiDcM++8CkSbDrrnDSSVFBP21aqaMquXwq4H8LjAKGA+cl/88tTFhN1HffxV0jO+4YY1a/8EKMza7xq0UaltVXh8ceg1tvhZdfjpO/W29t0g0d86mA/wOwJfCBu+8IbAHMLEhUTVFlZVw2X3UVHHdc3EXSp0+poxKRmphFA8eJE+O3e9RRUSz96aeljqwk8kkm8919PoCZLefubwEbFCasJmTBgqjQ69077hYZPjwaJLZsWerIRCQXXbpEf15DhsRdXz16xMBbTUw+yWSGmbUBHgX+a2aPAR8XJqwmYvJk2GabuNXwgAOiW5Rddy11VCKSr2bN4I9/jCKvLl3iVv5DDoGvvip1ZEWTUzIxMwNOdPev3f1c4C9Ej8F7FTC2xuv99+OSeJNN4vHQoXD33RqXXaTcde8ejYrPOQfuuw/WWSdGd/ym8fc8lVMycXcnrkiqnj/n7sPc/YeCRdYYTZsGgwZFh3L33BN1I5MmxVmMiDQOyywTpQ0TJsAOO0QX9127Rn9fc+aUOrqCyaeYa6yZbVmwSBqzDz+MFusujCkcAAAP4klEQVTrrRf3qR9zDLz7Llx9ddwVIiKNz2abxR1flZVRnH3mmZFULrkEvv221NHVu3ySyY7Ai2b2btIC/nUzm1iowBqF6dPh2GOhWze4/XY4+uhIItdcAx07ljo6ESmGXr1i8K1x42CrraITya5d4bLLYjyiRiKfZLIbsC6wE9ECfkDyX7J99FF0DtetG9xyS9SPTJ0K110HnTqVOjoRKYWttorOWV98MRLMqadGUrniCpg3r9TRpZZPMvm9u3+Q+Qf8vlCBlaWPP4YTT4y+s268EQ4/HN55B/7xD+jcudTRiUhD0Ls3/Oc/MGZMjEl08slRUT9kSFl3c69eg+vDp5/GbYHrrhuJ45BDIonceCOsvXapoxORhqhPn2iX8sIL0TblpJMiqVx9dVkmlbr2Gvy6eg0GPvssxl7v2jUaGg4cCFOmwM03x73mIiK12XZbeOYZGDUKNtpoyYnptdfC/Pmlji5n5rX0JWNmrYG2wMXAaRmz5rj7lwWMrVYVFRVeWVlZvDd0h/feg//+N778p56KL/uQQ+Css6KOREQkjZEjo53KqFHQoQMMGAD9+sHOO0O7dkUPx8wmuHtFrcvVlkwyVrg8sDfQBWhRNd3dz69jjKkVJZnMmhVdJVQlkKreQddaC3bbLa5M1luvsDGISNPiHknl2mvh2Wdh9uyYvsUWMbZKv36w3XawwgoFD6UQyeQ/wGxgArCoarq7X1HXINMqSDL57rsow3zmmfh75ZX4Ylu3jq6m+/WLL7NbN410KCKFt3BhNIB85pk4qR0zJvr0W265SChVx6TNNy/IeEeFSCZvuHuP1JHVo3pJJosWwauvLrnyeOEF+P77aMW6zTZLzgJ69YIWLWpfn4hIIc2dC88/v+SYNTFp7rfKKnHCW3XMWmedenm7XJNJPkfHMWa2ibuXf6X799/HwFPPPBOXkF8mVT+bbhpdnPTrF90gqOdeEWloWraE/v3jD+JGoGefXXLl8tBDMX2ddeJYNmAA/KrwTQLzuTJ5E+gGvA98DxjRbdemOb6+P3A10Bz4p7tfUsNy+wJDgS3dfamXHXW+Mlm0KCqyWraMLL7LLpHRO3TIf10iIg2FO7z99pKrlhEjYrC9Rx+t/bU1KMSVSZ3blJhZc+A6oq3KDGC8mQ1z9zezlmsFnAiMq+t75aR5c3jzzUgeqvcQkcbCDDbYIP6OPz7qW74szk23OSeTpMV7XW0FTHX39wDM7H5gT+DNrOUuAC4DTk7xXrlZY42Cv4WISEm1aAGrrVaUt8pnDHgzs4PN7Ozk+VpmtlWOL+8ITM94PiOZlrn+LYDO7v5ErjGJiEjDkE93KtcDfYADk+dziKKrXFRXlvRjZY2ZNQOGAINrXZHZ0WZWaWaVM2dqCHoRkYYgn2SytbsfB8wHcPevgGVzfO0MILOnw078dMjfVkAPYKSZTQN6A8PM7GeVPu5+k7tXuHtF+/bt8whfREQKJZ9ksiCpSHcAM2sPLM7xteOB9cysq5ktCxwADKua6e6z3b2du3dx9y7AWGCP2u7mEhGRhiGfZPJ34BFgNTO7CHiB6K+rVu6+EDgeGA5MBh5090lmdr6Z7ZFnzCIi0sDk3M4EwMw2BHYm6kCedffJhQosF0Xv6FFEpInJtZ1JPndz3QF86u7Xufu1wKdmdmuaIEVEpHHIp5hrU3f/uupJUgG/Rf2HJCIi5SafZNLMzNpWPTGzVcivBb2IiDRS+SSDK4AXzWxo8nw/4K/1H5KIiJSbfLpTudPMKoGdkkl7Z/etJSIiTVPOycTMlgM2B1ZOXrevmZV0pEUREWkY8inmeowlIy1+X5hwRESkHOWTTDq5e/+CRSIiImUrn7u5xpjZJgWLREREylY+VybbAYebWZ1GWhQRkcarKCMtiohI41askRZFRKQRy6sFu5ltBmyfPH3e3V+r/5BERKTc5NPR4x+Ae4DVkr+7zeyEQgUmIiLlI58rk6OI0RbnApjZpcCLwDWFCExERMpHPrcGG7Ao4/kiqh/bXUREmph8rkxuA8aZ2SPJ872AW+o/JBERKTe1JhMz6was7u5XmtlIor2JAScCHxU2PBERKQe5FHNdBcwBcPeX3f3v7n41MC+ZJyIiTVwuyaSLu0/MnujulUCXeo9IRETKTi7JZPmlzFuhvgIREZHylUsyGW9mg7InmtlRRHf0IiLSxOVyN9cfgUfM7CCWJI8KYFng14UKTEREyketycTdPwO2MbMdgR7J5H+7+/8KGpmIiJSNfDp6HAGMKGAsIiJSpvJpAS8iIlItJRMREUlNyURERFJTMhERkdSUTEREJDUlExERSU3JREREUlMyERGR1JRMREQkNSUTERFJTclERERSUzIREZHUlExERCS1oiUTM+tvZlPMbKqZnVbN/JPM7E0zm2hmz5rZ2sWKTURE0ilKMjGz5sB1wG5Ad+BAM+uetdgrQIW7bwo8BFxWjNhERCS9Yl2ZbAVMdff33P0H4H5gz8wF3H2Eu89Lno4FOhUpNhERSalYyaQjMD3j+YxkWk2OAp6qboaZHW1mlWZWOXPmzHoMUURE6qpYycSqmebVLmh2MDHG/OXVzXf3m9y9wt0r2rdvX48hiohIXeU8bG9KM4DOGc87AR9nL2Rm/YAzgV+4+/dFik1ERFIq1pXJeGA9M+tqZssCBwDDMhcwsy2AG4E93P3zIsUlIiL1oCjJxN0XAscDw4HJwIPuPsnMzjezPZLFLgdWAoaa2atmNqyG1YmISANTrGIu3P1J4MmsaWdnPO5XrFhERKR+qQW8iIikpmQiIiKpKZmIiEhqSiYiIpKakomIiKSmZCIiIqkpmYiISGpKJiIikpqSiYiIpKZkIiIiqSmZiIhIakomIiKSmpKJiIikpmQiIiKpKZmIiEhqSiYiIpKakomIiKSmZCIiIqkpmYiISGpKJiIikpqSiYiIpKZkIiIiqSmZiIhIakomIiKSmpKJiIikpmQiIiKpKZmIiEhqSiYiIpKakomIiKSmZCIiIqkpmYiISGpKJiIikpqSiYiIpKZkIiIiqSmZiIhIakomIiKSmpKJiIikVrRkYmb9zWyKmU01s9Oqmb+cmT2QzB9nZl2KFZuIiKRTlGRiZs2B64DdgO7AgWbWPWuxo4Cv3L0bMAS4tBixiYhIesW6MtkKmOru77n7D8D9wJ5Zy+wJ3JE8fgjY2cysSPGJiEgKLYr0Ph2B6RnPZwBb17SMuy80s9nAqsCszIXM7Gjg6OTpt2Y2JUVc7bLX34Q01W1vqtsN2nZte92snctCxUom1V1heB2Wwd1vAm6ql6DMKt29oj7WVW6a6rY31e0Gbbu2vbCKVcw1A+ic8bwT8HFNy5hZC6A18GVRohMRkVSKlUzGA+uZWVczWxY4ABiWtcww4LDk8b7A/9z9Z1cmIiLS8BSlmCupAzkeGA40B25190lmdj5Q6e7DgFuAu8xsKnFFckARQquX4rIy1VS3valuN2jbm6qibLvp5F9ERNJSC3gREUlNyURERFJrEsmkqXblksN2n2Rmb5rZRDN71sxyup+8HNS27RnL7WtmbmaN5rbRXLbdzPZPvvtJZnZvsWMslBz2+bXMbISZvZLs97uXIs76Zma3mtnnZvZGDfPNzP6efC4TzaxnvQfh7o36j6jwfxdYB1gWeA3onrXM74EbkscHAA+UOu4ibfeOwIrJ42Mbw3bnuu3Jcq2AUcBYoKLUcRfxe18PeAVomzxfrdRxF3HbbwKOTR53B6aVOu562vYdgJ7AGzXM3x14imjP1xsYV98xNIUrk6balUut2+3uI9x9XvJ0LNH+pzHI5TsHuAC4DJhfzOAKLJdtHwRc5+5fAbj750WOsVBy2XYHVk4et+bn7d3KkruPYunt8vYE7vQwFmhjZmvUZwxNIZlU15VLx5qWcfeFQFVXLuUsl+3OdBRx5tIY1LrtZrYF0NndnyhmYEWQy/e+PrC+mY02s7Fm1r9o0RVWLtt+LnCwmc0AngROKE5oJZfv8SBvxepOpZTqrSuXMpPzNpnZwUAF8IuCRlQ8S912M2tG9Ex9eLECKqJcvvcWRFFXX+Jq9Hkz6+HuXxc4tkLLZdsPBG539yvMrA/Rtq2Huy8ufHglVfBjXFO4MmmqXbnkst2YWT/gTGAPd/++SLEVWm3b3groAYw0s2lEGfKwRlIJn+v+/pi7L3D394EpRHIpd7ls+1HAgwDu/iKwPNERYmOX0/EgjaaQTJpqVy61bndS1HMjkUgaS7k51LLt7j7b3du5exd370LUF+3h7pWlCbde5bK/P0rcfIGZtSOKvd4rapSFkcu2fwjsDGBmGxHJZGZRoyyNYcChyV1dvYHZ7v5Jfb5Boy/m8obblUtB5bjdlwMrAUOT+w0+dPc9ShZ0Pclx2xulHLd9OLCrmb0JLAL+7O5flC7q+pHjtg8GbjazPxHFPIc3ghNHzOw+otiyXVIfdA6wDIC730DUD+0OTAXmAUfUewyN4HMUEZESawrFXCIiUmBKJiIikpqSiYiIpKZkIiIiqSmZiIhIakomIiKSmpKJiIikpmQiApjZIjN7NRnf47VkrJdmybwxS3ldGzP7ffEiFWmY1GhRBDCzb919peTxasC9wGh3P6eW13UBnnD3HgUPUqQB05WJSJakn7KjgeOTvoy+BTCzlmb27+TK5Q0z+w1wCbBuclVzebLco2Y2IbnKOTqZ1sXMJpvZzcn0p81shar3NLNDkxHwXjOzu5JpB5vZS8m6bzSz5tmxJqMG7pI8vtDM/l7oz0ekOroyEeGnVyYZ074CNgTedfeVzGwfoL+7D0rmtwbaknVlYmaruPuXSbIYT3Tt34roF6nC3V81sweBYe5+t5ltDDwMbOvus8xsFWB1YuCuvd19gZldD4x19zuzYtwBOB+4GRhIdFi5qN4/IJFa6MpEpGbZY0C8DvQzs0vNbHt3n13D6040s9eI3og7s6R79/fd/dXk8QSgS/J4J+Ahd58F4O5fEj3b9gLGm9mryfN1st8oGWHPgJOAA5RIpFQafa/BInVhZusQPer+2DW/u79tZr2I3lcvNrOngewrhb5AP6CPu88zs5FEN+cAmePFLAKqirmM6gdsu8PdT68lzk2ANYBZ7j4n5w0UqWe6MhHJYmbtgRuAazO7JzezNYF57n438DegJzCHKMKq0hr4KkkkGxIDb9XmWWB/M1s1eZ9Vkmn7JjcDYGarmNnaWXGuAdxDjO8918x+WacNFqkHujIRCSskxUnLAAuBu4Ars5bZBLjczBYDC4Bj3f0Li7HU3wCeAs4CjjGzicQIhmNre+NkzI2LgOfMbBHwirsfbmZnAU8ntygvAI4DPgAwsxWJepbB7j7ZzC4ALiXG8hApOlXAi4hIairmEhGR1JRMREQkNSUTERFJTclERERSUzIREZHUlExERCQ1JRMREUnt/wE5PhIaBIjtZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "makePlot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Another solution scheme could preserve the values of the concentration.  I choose not to do that here for simplicity.  Doing so would permit a small animation of the diffusion process, however.  This is left to the student to implement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Reading Assignments and Practice\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Analyze the spectral difference schemes in the following two examples.  For reference and attribution the codes can be found [here](https://open.umich.edu/find/open-educational-resources/literature-science-arts/parallel-spectral-numerical-methods).  The materials are provided under a Creative Commons license with attribution to the original authors whose names can be found at the above link."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# %load Heat_Eq_1D_Spectral_BE.py\n",
    "#!/usr/bin/env python\n",
    "\"\"\"\n",
    "Solving Heat Equation using pseudospectral methods with Backwards Euler:\n",
    "u_t= \\alpha*u_xx\n",
    "BC = u(0)=0 and u(2*pi)=0 (Periodic)\n",
    "IC=sin(x)\n",
    "\"\"\"\n",
    "\n",
    "import math\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator\n",
    "\n",
    "# Grid\n",
    "N = 64; h = 2*math.pi/N; x = [h*i for i in xrange(1,N+1)]\n",
    "\n",
    "# Initial conditions\n",
    "v = [math.sin(y) for y in x]\n",
    "alpha = 0.5\n",
    "t = 0     \n",
    "dt = .001 #Timestep size\n",
    "\n",
    "# (ik)^2 Vector\n",
    "I = complex(0,1)\n",
    "k = numpy.array([I*n for n in range(0,N/2) + [0] + range(-N/2+1,0)])\n",
    "k2=k**2;\n",
    "\n",
    "# Setting up Plot\n",
    "tmax = 5.0; tplot = 0.1\n",
    "plotgap= int(round(tplot/dt))\n",
    "nplots = int(round(tmax/tplot))\n",
    "data = numpy.zeros((nplots+1,N))\n",
    "data[0,:] = v\n",
    "tdata = [t]\n",
    "\n",
    "for i in xrange(nplots):\n",
    "    v_hat = numpy.fft.fft(v)  # convert to fourier space\n",
    "    for n in xrange(plotgap):\n",
    "        v_hat = v_hat / (1-dt*alpha*k2) # backward Euler timestepping\n",
    "\n",
    "    v = numpy.fft.ifft(v_hat)   # convert back to real space\n",
    "    data[i+1,:] = numpy.real(v)   # records data\n",
    "\n",
    "    t = t+plotgap*dt    # records real time\n",
    "    tdata.append(t)\n",
    "\n",
    "# Plot using mesh\n",
    "xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,\n",
    "        linewidth=0, antialiased=False)\n",
    "fig.colorbar(surf, shrink=0.5, aspect=5)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# %load Heat_Eq_1D_Spectral_FE.py\n",
    "#!/usr/bin/env python\n",
    "\"\"\"\n",
    "Solving Heat Equation using pseudo-spectral and Forward Euler\n",
    "u_t= \\alpha*u_xx\n",
    "BC= u(0)=0, u(2*pi)=0\n",
    "IC=sin(x)\n",
    "\"\"\"\n",
    "\n",
    "import math\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator\n",
    "\n",
    "# Grid\n",
    "N = 64                     # Number of steps\n",
    "h = 2*math.pi/N                 # step size\n",
    "x = h*numpy.arange(0,N)    # discretize x-direction\n",
    "\n",
    "alpha = 0.5                # Thermal Diffusivity constant\n",
    "t = 0\n",
    "dt = .001\n",
    "\n",
    "# Initial conditions \n",
    "v = numpy.sin(x)\n",
    "I = complex(0,1)\n",
    "k = numpy.array([I*y for y in range(0,N/2) + [0] + range(-N/2+1,0)])\n",
    "k2=k**2;\n",
    "\n",
    "# Setting up Plot\n",
    "tmax = 5; tplot = .1;\n",
    "plotgap = int(round(tplot/dt))\n",
    "nplots  = int(round(tmax/tplot))\n",
    "\n",
    "data = numpy.zeros((nplots+1,N))\n",
    "data[0,:] = v\n",
    "tdata = [t]\n",
    "\n",
    "for i in xrange(nplots):\n",
    "    v_hat = numpy.fft.fft(v)\n",
    "\n",
    "    for n in xrange(plotgap):\n",
    "        v_hat = v_hat+dt*alpha*k2*v_hat # FE timestepping\n",
    "\n",
    "    v = numpy.real(numpy.fft.ifft(v_hat))   # back to real space\n",
    "    data[i+1,:] = v\n",
    "\n",
    "    # real time vector\n",
    "    t = t+plotgap*dt\n",
    "    tdata.append(t)\n",
    "\n",
    "# Plot using mesh\n",
    "xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,\n",
    "        linewidth=0, antialiased=False)\n",
    "fig.colorbar(surf, shrink=0.5, aspect=5)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('t')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}