Fourier-Spectral.ipynb

{
 "metadata": {
  "name": "Fourier-Spectral"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Fourier spectral methods in Matlab (and Python)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.\n",
      "\n",
      "See <http://faculty.washington.edu/rjl/classes/am590a2013/codes.html> for more IPython Notebook examples."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These examples are based on material in Nick Trefethen's book Spectral Methods in Matlab.  The m-files for this book are available at <http://people.maths.ox.ac.uk/trefethen/spectral.html>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymatbridge\n",
      "ip = get_ipython()\n",
      "pymatbridge.load_ipython_extension(ip)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Starting MATLAB on http://localhost:54980\n",
        " visit http://localhost:54980/exit.m to shut down same\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "MATLAB started and connected!\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you want to run this notebook you will have to set this path properly..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "path(path,'/Users/rjl/SMM/all_M_files')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Program 5"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This example is directly from p5.m found at  <http://people.maths.ox.ac.uk/trefethen/spectral.html>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "\n",
      "% p5.m - repetition of p4.m via FFT\n",
      "%        For complex v, delete \"real\" commands.\n",
      "\n",
      "% Differentiation of a hat function:\n",
      "  N = 24; h = 2*pi/N; x = h*(1:N)';\n",
      "  v = max(0,1-abs(x-pi)/2); v_hat = fft(v);\n",
      "  w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* v_hat;\n",
      "  w = real(ifft(w_hat)); clf\n",
      "  subplot(3,2,1), plot(x,v,'.-','markersize',13)\n",
      "  axis([0 2*pi -.5 1.5]), grid on, title('function')\n",
      "  subplot(3,2,2), plot(x,w,'.-','markersize',13)\n",
      "  axis([0 2*pi -1 1]), grid on, title('spectral derivative')\n",
      "\n",
      "% Differentiation of exp(sin(x)):\n",
      "  v = exp(sin(x)); vprime = cos(x).*v;\n",
      "  v_hat = fft(v);\n",
      "  w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* v_hat;\n",
      "  w = real(ifft(w_hat));\n",
      "  subplot(3,2,3), plot(x,v,'.-','markersize',13)\n",
      "  axis([0 2*pi 0 3]), grid on\n",
      "  subplot(3,2,4), plot(x,w,'.-','markersize',13)\n",
      "  axis([0 2*pi -2 2]), grid on\n",
      "  error = norm(w-vprime,inf);\n",
      "  text(2.2,1.4,['max error = ' num2str(error)])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAbBUlEQVR4nO3d2ZKj\nyLIFULh2/v+XuQ9U0RSTkJhiB2tZW1uWpFS6CAcnBkHbdV0DAGn+7+kAAOAXChgAkRQwACIpYABE\nUsAAiKSAARBJAQMgkgIGQCQFDIBIChgAkRQwACIpYABEUsAAiKSAARBJAQMgkgIGQCQFDIBIClht\n2rY98os//zo8ZU/Srr3m4+/u3COGl9mD7qSA8Y+u654OAQry7R5hD7pTa3PXZDj767qubf80bv/D\n+MRweHz84o1nm3/PK+UMJ5qk1sdUbFZSd/h57U3Gebv2Ryd/Yvyew+Pz91+MYfH19p1z/e/pADjT\nuG4tPtuMhgrn+9iwp82fbZYqGZxikrGLhar5N1GHF0weGedqs5Kuk3cYl5z5WzWj3WqyLzT/1rm1\n16+9P8cpYDgr5El7Oihrp03fvn7DV79y9fuzkwKGs0KeNOkPzV+w8exi6n58w6/+xA8BH3x/9lPA\naraxt2yM+689C1eYTCk16/Ovw/DdYuquTXF9/It7sn0+dDl/zeIo4s735zdOvYGCGA9gP8voAYjk\nZAeASHpgAERKLWBrk6iD+0OCB8l5XihvFeLH9UW3RQIlULp4rbweWNd1G1VK94u32d4joGJ5PbBt\nk+u1KGYcl1geZD7HlZ/5VRWwxc1dfhssCv02TH1h51aCyhqicKGRp2d+3hDimojNDcBZauiB9ScR\nrtcC8CqR3d79Qvv1lCM0hULDphwRKVTPECIAr6KAARBJAStU6JoUYXNEbkOERh4a9kABAyCSAlao\n8qdPFwmbI3IbIjTy0LAHChgAkRQwACIpYIUKnVwVNkec3hBt29zTtqEpFBr2QAErVOjYtLA54tyG\nGA7ONxylQ1MoNOyBAgZUKLxrwS4KGFCh8K4FuyhghQodmxY2R1zUEDcUs9AUCg17oIAB1dIPq1sN\nt1OpUujkqrA5Yrshxr2Fjy3WtrdWr9AUCg17oAcGPKNf475zEOvnsa6us6CjWhX2wCJuYwPfWrtf\n63gaIzTzT+8t3dz94ilV9cDatk2fkxzc80FO/5pn6PYvP+zhtuPjO48PhqceiW1s/4ZcfOUpDfFI\n9So/hRaFhj2oqgfW78DpTdK74WA0/prnWX+thGPoD0LDHvQ5//in6DOq//9vsax9hEmKFriLP77x\nfxMa9qCqHtiicbdsXNsm3bUXPjXWttNXlhDhs08FdegXe2Y3Z/5aj6p/6u9///3cNOORz38GA2Z/\n67/q1T81TGsthtG/+N+A/vlXsSlXyFNBmV/hdFE7mgNrzYet+2pN12sVkkJ7svrZzN95xNsOqi8/\nw1v1P6/9yuJTa683K/atQjJ/W/09sFBXnwFdNCaTcuI2ERp2U1Lk43Tquv/++9Kkq7T1DvPlhQ9W\nqXIa4iuhYQ+qmgNjp2E/n9ew4k+5Xmo8QjjM9U5GDks4X56HcOkq9qGLpoP1TgGdxCMiesH3GM0N\nbO3qO1/2HqEp9MgQ4nZeNbtPj74d3N75ein9lYjM1wOr32huvGk+HRHG58t2eM7yVSINSfht9eJt\nzIGxxdGBR/w0f8brKGCFumhy9duDwrevD50TDg27HGd11n9oiELqXGgKhYY9MIRYqLNGnw/m57cH\npvIHzReFhl2f3xpi5y9dutYjNIVCwx4oYDX7bXed/IrViWwzV8pTDCHW5vQlGP1sxGQlCMDjFLBC\n/TA2PVyMp//h3JPi8cV+NkILHVIPDbsE52ZabkOERh4a9kAB4xfhaQ/UQAEr1MHJ1afmJELnhEPD\nrs/VDXHdZUFCUyg07IECVokbukThqc75LN/gWVYh1uC248jk2okOXsCD9MAKtX9y9alCMrmH099g\nIifHQsN+1hWJl9sQoZGHhj3QAyvUxtj0JOUe7AYNl/0dPRLZKQsN+ylXXl3+8oa46OvMoSkUGvZA\nDyxbaekXfj7HZ5P7dcGDFDAAIilghUofm85iaxcityHKifyrQMoJ+zd5c2Ab958dN0b62O6iMi9L\n+O/SxIjb4FXl+B2Zh/3m4xvUcb/TWu/g/Ntl5KI3RVgBGx8dF4+U1Rw7l8pzQJ79vdt9//PDweyX\nmzYf94h/X/znh8nXIcY/r73BpE2jr+k+viLaWX+w8BRab/o/O+xSSnTlH3PCCthH/dlo4cm0Uzmr\nDb+1tDrxqVj444fFF1W2YPiY2V4fP+bsCzCrLyu56WubA+u6ruu68Vhi27bDP9ceL/Opfz9XiRFu\nPzWO/9+jZykRbj81eUGctczv9TcZmB+b+u/2LX7Dr5ymOfLU7G5BxUX421OzFm6bpp038dC4k8f/\nvrJtmqScD5uvaPcNmAxPbbymcG3b9r37QcTn+LeB/nmq5Pj35FKZ1vaIxbDX5ro+Di3Of+Ui92zt\nK3qWz+bJdj9p7fN++3iB6hlCLPxA862u64LSaBDaBKFhf2vtUy4+ft11bzfc0xBXrEO5OYWGCbw9\nM3lrL9vOh4gjalgBGw8PTk42F58KFbcIYtH8ND/9ExXo0rTXXmUaTiz271O1NmVYAWuW9tLhkfS6\n1av1QL/zbJFv1ZH28IPaFnEkGk+bt/99zyZpKnWwHXY/Uby2TOBBoVu7Pnc2xLljpKEpFBr2QAEr\ny9tOpsN3H+BJeUOILxE6LiRsjshtiJsjf8n3rz/SAwMgkgL2pGHSqz8NCj8Z2mv8eR9Zqw096ZdO\nAXvMuHQ1SxcIuD+k4/aHPf7gj3/W0K1dn9yGuC3yc1cp527wnjmw+3x1aYPQsenfwp7XsJs/fejW\nrk9uQ4RGHhr2QA/sGeFpc761ixsBg1q/JPozBQx4rxJGsPmZAlao0LFpYXNEbkPcEPkV3a/cDd5T\nwO4wuX+BQYBFVicCX7GI43K/la7QydXjYU9WJ96zGUK3dn3KbIg994W4OvKL9oUyN/h+Ctj5gu6D\nVbjHVyfyBneeKnEuQ4jXslccZHUiiYq6VnXFFLBChU6uCpsjchtiHPn4fl27f33rxdd1EHM3eK+q\nIcQSbmh5Vj6Ejk3fcz/40/9IxNZeS+/xMSjig2woNv4hrrX0OxL5DzeoPEuxG3ynegrY+AbY2zfD\nXpuSPfj4+G6N7tx4ovkN4Cf3U2/WRxrPeLwrYYJkO73TD0OP25gGmzx+54RZeO/oDoYQzzG/sKFD\nyrkmB5HJ+o47fy5Q27bpY0HVmKw6/jgZNnm2wNu9lqz+Ajbet+e7+Vq6fPt4100PIpOfv32q/eus\nN7znqUn894Tx95Hl3f6sx0vWdV3XdZNtMmmOxcdLe2r4fzERjvv9wxBunyELbzKEPzqFbUd3IW/H\nqdX+ezOK5u+ZWdf9+a1F537kxV148mDJtobasrRLYyzt0ljiuGlOHEI812Lk5bsz7Pk5xPbjTWFN\nvGFy+OiL08cR8j2vKdyDYbdLA4OLDy4+1e6YNdjIzI+vv24Rx9oGj0iheubA9pvMV53++ElBlp46\ni+4Mez43tv14sz6B8e3jV9u/GSOOMvs9+Fnmbb3d9JNB7G8D3/P6jUw+S3ry1FPAxkMoH1tl7fmz\nHuc21zdl2zTPN/M8vfvS9VXak0irbqingDV2YKq2tvJQ2l/hin7PpWM271T/Io5QKZOoE8LmiEIa\n4ofqtTPy0tYnF7LBf1ZVD6wmoafVwuaIZxvi2ynPSZf4/ICul575emAATfN31V94n+RdFDAAIilg\nhQodmxY2RzzYEAf/cmgKhYY9UMAAylpbwU4KWKFCJ1eFzRGPL+Jofl0oGJpCoWEPFDCAP8KP56+j\ngAEQSQErVOjkqrA5IrchQiMPDXuggBUqdGxa2ByR2xChkYeGPVDAAIikgAEQSQErVOjYtLA5Irch\nQiMPDXuggAEQKe9q9Bu37xufTaRPTobGL+yn1HFr5tyPEBp5aNiDsB7YcBfa8Y1ox4Zn74/tXKFd\ne2Hfr23b6PjHcj9IaOShYQ/CCthHNe3MsEcdZ2zwg7whxG39njweTsmtZ6GRC7sQoZ8oNOwmNvLQ\nsHtFF7DJlv14mjl/gTNTarJ/j5D5vEHRBeyrnbCOSWzYIMNhrOgCNjdeuzEeJJws67CfA1RPrwWA\nSLWtQgTgJcKGEL+SO6KYG3kTOBkZvbUXhX6i0LB7cWnfhG/wXrUFbLKSPq6F5t8HiBC3JLcPOHRr\nL4rO/NCGiEv7ppbMN4RYotBkyt0N+u+/hwZfjdDtH5056ZlfbQ+sAtGJFSS3v1IlrXCbCjJfD6xE\niadF/YjE+P/wFWnPt/TACpW1GzdVnM3xuLjMkfbPqnmjh66x+fYCWqWJ25ND82RD4ieS9vdLzJOJ\nvI0OAI05MABCKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQED\nIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQ6X9PB3Cmtm37H7quezYSuI2057Vq64F1Xdd1\n3bBLwxtIe96pqh7Y/AzULs1xhfdsFsOT+RxXeOY3lRWwZmm/Lb8NFrVtmxh5fWGnVIL5R6isIQoX\nGnl65tdWwPrGiNj0cIo+2xOPnnBQPXNgihavpXrxTvX0wMaT2PZnXqLPeZnPO9VTwJq69t7QzyLs\nm+VGvij344RGHhr2oJ4hRABeRQG7Vds2puoATqGA3WcoXXtqWOiaFGFzRG5DhEYeGvZAAStU6Ni0\nsDkityFCIw8Ne1DVIo5y9Kc1Q27Mz3LWXhCeTgD3UcDONx8qXCtU4d134I+1bzKMx+gml1lYvOrC\nka8DjS+r8ZJvVihg15okz+I/F8tYfVemKVlo2GdZGwC4f2BgrSHGh+MCD83ji6FsX9Nr/srJs5Mi\ntPMzLk5lrcUz+cVyNuMPFDB4tfGhr22Xx73Hjz9lcjiefIN78oL5cXmx3zOUw42+0fHj+/x9Ft95\niHztHTaCmXfmosvSfgrYyfpdfTLFtW38+uFIEZp/wk737LD2tw0xHzFbtNhDGn533lOZv2b+Jxa7\nWXt6Qjs7WMNbbXfvPn7qbemZr4CdaVR+vvvF4fXjGgZXmxzi1oYQi7V494mP45C/2TjWj3ty23Vu\nzbi8Lb7VpJ7tjGr++sooYKc5q/CoYVxkPKc1Xwc7T7nxI18NKtxjcc3C2hDizr7a9p9bfMNmvWP0\nwwzT4q98OyuWPrO1nwJ2jnNLTm4NC91zQsP+yvac1sdPPyw4mvxwrm8bYlJU1n59z9KPj6/Z0wNr\nZsVm8vjGK9fe6od1K/t/JT3zs6P/6NLmGe871+zM174/e4Tu4UurGP55wc+f6az3oXARme9KHABE\nUsAyVD0Ryx267r8VRsWfWMMuClihZjPGTwXyndAlT6Fhj+38CMcTafwOp6dlbkOERh4a9sAijkMu\nrSvDbDlsuP8YtL12kXTD+o/y21cP7EdXN+1oVdJ/i57LV/6s76LQsHuTyxHd+6dPzszchgiN/GPY\nhR95FLBf3H9iElTDKER/Yb3h58XHPeWpjafKF7BQ8oiLVoI+1bMuv0dfn4jFxM0szjvDlpaVCfoC\njx7Y1+7ZXbPOgwbCfqEThwdyGyI08knYbfvn+NZfharw6tUoYN+67WRz5VoApQ8kRnRW5kLDrk9u\nQ4RG/m+vPe8rFlYhfqGEoZLcq0xxqWcPoNIyUQUXVamqgA3d4XN35gI7PeOrskIJxsMD0jJOaJPV\nNoTYLykOHY8eC/0IwqY5cM6X2xC5kUerqoCt3dHg4GLitXf7+Q2PPzWPqrQIK3gqbknxs2yqLHW0\nV8YS4a+0S3dWPfaG//1cwtYqLZ66pSyjn7g/7AomVN5jz4RlROZXOAd2xUYvqh3HF/KBQpgDK9b8\nQFFNA1VVwJoLqlfJa6ss/aIorpFYvsrapZ45sL77de7UxYP74c6PUNo3w0InjULDrk9uQ6REPgkz\nJew19fTAyh+u/UroxxE2zYGxgdyGKDby7bnJYsPeqZ4CdrqUYRADiUAzW97V/7Pucd16hhDPldXY\npQ0kgpx81vy6UEEHtP0UsEKFjk0LmyNyG6K0yHeWq9LC/pYC9o++NbO6X73hhDc8IamHTth1+svG\nr/3zPcyB/aeoAvDb5OrjBTh0Tjg07PrkNsSdkQ/HqPkX7ybzXh/lbvCeHlidCinDvJxO2NXmE13h\nJek7emB1elUSQ5UWVxUypgdWqNDJVWHfr/DrDn915C32U3x0deTDqsJxH+v4eWruBu/pgf1R2sKN\nH8amx7/x1McJHVIPDbu54NLVF9k5N1Ns/B8djHy+fT59AfnIXxu/T+oG7ylgwE1KO00sxHxRRmPM\ncB9DiE1T434l+xmPKz56Z7XP4T0d4a1bY1SuFjbN3xHC9t+f3QlvWbljDqfYM6hSZvU6ZTjo/o9W\n8ijWho2wC/9Ea0OIRYW9/w52RYX9lZ2HmsHaKWY5O2xEW+iB1Uw/jBKcuOigGhctynibt8+Bldn9\namInV4V9s67rhtGekj9FH9rH3a3kj7BtHPm8x/nxlpKuPPCblxaw99w6drhWfbGlmoPSj0F1Gy8v\nNBxyujcWsBem0eOXmILmlbf+qf568M963RxYSvW6aBXQ1R8/ZfHSRGjY9YlriGFVYR953EV14zb4\nxOsKWMpJ0EXjQld//NDhrNCwE20vLMpqiH+/v9WN12X0yv80WRt87o1DiAAH3XOlDLa9sYC9LbdK\nuMQUDF44E8ZFXjeEmOKisemrvxkWOqQeGnZ9im2If6+gMS3AXVdu5NtCwx68sQf2cs5/eVxWEk6u\nVRg0xVU9PbBCXTq5el0/LHROODTs+hTYENtzXaPHi4t8j9CwBxUWsPRO8T1cZYpnpWRg+BG+clUN\nISpd33rPFUl+MyyPbmyiaxS+bR1RCldVD6zruvQe8eDmYnzWXws9h1gMO/OjxBhv3u07jNxpEtX4\nervjL3gt/WJkuoSGPaiqgC0a39um1FsELTzVX6T10r81ccrfGt/O4/FtuP+p8SVx57/FbZ49AR2u\nuDZfZ/hR6KlzaNiDgDu+fKst9a5Ipdk5O/1a0dsnIvP33yfsBtHNfYWIFKpqDoyvTJIzaFnzDYat\nEbEbhxrfKuFxhYTBV+ofQgx1/xDWKTtw6MjbJGy1/Dazs6jH8ufgX64j8+NU2ANzvvyzrK+XXiR0\nC0T3FB9PvNBGp8ICVoenDkaTfti3UcQeQ7vcbxSkn0T3/tawm7b+aIj4z18/KDTzQ8MeKGBsecOZ\nafrxvz8G1VHG7jGsNmwCz1cYMwcGdYr7AsnQ+7/0b00KfbFb48Gngr5GEjxuvkfuxMCDkR9ZT5y4\nwbc/b1GfaHJYGQc2ibOosPdr2/aGYa0rlu/nbvC1sCM+kSHEQj2YOvPl9fMH13+39IwfixtEytq8\nP+i67upT/4veP7RpQsMeKGB8MMx1j9d3hKa9L6uW7+Ayog1x5yt8ZA6MXSaHlZAR8n/UXb3ST6Uv\n0v69EWX798KG1EQBK1TKJOqEsDni3NUDw1vdsJg2NIVCwx4oYPxotprroThWTPqLpYXHmvNWVZzz\nPpTMHFihChwRGk+G9eY1rISwx+fdvXFQiyfjJYRN87chxhfmSJm4Ck2h0LAHemB8Z5zw80mFAk97\nJxGG77BvMTRTn2Pzk5KJcSd7fjMUjV4rPTBONj5lvnrV4uT9CyyfnGJewya96sWn1K3q6YEVKmhy\n9d/DRDucMm+sWmzb5XmptcmqxdfPZ7nmt83deQgL2tp1298QQ0oUsrg0NIVCwx4EfNf6iIgvk9dq\n566xOKO2eH69833OFZpCoWFv25Mh1X3ox0SkkCFErrL2xec9yxc3itbiWxW/o3EHafA2ChgX+nhA\n2Shsa2sdv3p/arLW3NLgtcyBFSp0bHpP2MM01cb6wMW1jtddSSF0a9cntyFCIw8Ne6AHxjMWS9FG\nfXKWDUzogRWq/OnTRcLmiNyGCI08NOyBAgZAJAUMgEgKWKFCJ1eFzRG5DREaeWjYAwWsUKFj08Lm\niNyGCI08NOxBVasQh7OJ9FaB/aQ9r1VPD6y/8EkvvV/cxHbthf0Iaf+40MhDwx7UU8DgnXS8eK2q\nhhAX5Z5ihEYu7EfML70a+olCw25iIw8Nu1d5AXNySk0mx5o+vfsHJ6ku83mDygsY1GStLClXvFPA\nHV/2sxyLt1nsk8FLVFXAAHiPmocQcztkuZE3ITdyHYve2otCP1Fo2L24tG/CN3iv2gI2zqfE3Brm\n57Mij1vRNF4BEbe1F0VnfmhDxKV9U0vm+x5YiUKTKXc3aNs2N/hqhG7/6MxJz/xqe2AViE6sILn9\nlSpphdtUkPl6YCVKPC3qRyTG/4evSHu+pQdWqKzduKnibI7HxWWOtH9WzRs9dI1N+jd74vbk0DzZ\nkPiJpP39EvNkIm+jA0BjDgyAUAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEU\nMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAk\nBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQ\nSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAA\nRFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoY\nAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJIC\nBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIik\ngAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAi\nKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyA\nSAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQED\nIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLA\nAIikgAEQSQEDIJICBkAkBQyASAoYAJEUMAAiKWAARFLAAIikgAEQSQEDIJICBkAkBQyASAoYAJEU\nMAAiKWAARPp/3YI65ES2xb0AAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Illustration of spectral differentiation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To make this a bit clearer, first illustrate how to compute the second derivative of periodic function.\n",
      "Start with $$u = \\exp(\\cos(x)),$$ and check that the numerical approximation agrees well with $$u''(x) = (\\sin^2(x) - \\cos(x))  \\exp(\\cos(x)).$$\n",
      "\n",
      "The only tricky thing here is the order of the indices in the wave number vector."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "N = 16;\n",
      "x = linspace(2*pi/N,2*pi,N);\n",
      "ik = 1i*[0:N/2 -N/2+1:-1];   % i * wave number vector (matlab ordering)\n",
      "ik2 = ik.*ik;                    % multiplication factor for second derivative\n",
      "\n",
      "u = exp(cos(x));\n",
      "u_hat = fft(u);\n",
      "v_hat = ik2 .* u_hat;\n",
      "v = real(ifft(v_hat));     % imaginary part should be at machine precision level\n",
      "\n",
      "error = v - (sin(x).^2 - cos(x)) .* exp(cos(x));\n",
      "norm(error,inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "text": [
        "ans =\n",
        "     3.909521448797193e-07\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Illustration of solving a periodic boundary value problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's solve the boundary value problem\n",
      "$$u''(x) = f(x)$$\n",
      "on $0 \\leq x \\leq 2\\pi$ with periodic boundary conditions and the constraint $\\int_0^{2\\pi} u(x) dx = 0$.\n",
      "\n",
      "Use $f(x) = (\\sin^2(x) - \\cos(x))  \\exp(\\cos(x))$ so the solution should be $u(x) = \\exp(\\cos(x)) + C$, where the constant is chosen so the integral constraint is satisfied.\n",
      "\n",
      "We now have to divide by ik2, with the complication that 1/0 should be replaced by 0.  This results in the $\\hat u_0 = 0$, which gives the integral constraint."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "\n",
      "N = 16;\n",
      "x = linspace(2*pi/N,2*pi,N);\n",
      "f = (sin(x).^2 - cos(x)) .* exp(cos(x));\n",
      "f_hat = fft(f);\n",
      "\n",
      "ik = 1i*[0:N/2 -N/2+1:-1];   % i * wave number vector (matlab ordering)\n",
      "ik2 = ik.*ik;         % multiplication factor for second derivative\n",
      "ii = find(ik ~= 0);   % indices where ik is nonzero\n",
      "ik2inverse = ik2;     % initialize zeros in same locations as in ik2\n",
      "ik2inverse(ii) = 1./ik2(ii);   % multiplier factor to solve u'' = f\n",
      "\n",
      "u_hat = ik2inverse .* f_hat;\n",
      "u = real(ifft(u_hat));   % imaginary parts should be roundoff level"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting the solution shows that it is a shifted version of $\\exp(\\cos(x))$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "plot(x,u,'b-o')\n",
      "hold on\n",
      "v = exp(cos(x));\n",
      "plot(x,v,'r-o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAATfklEQVR4nO3d6ZKj\nVhaFUejw+78y/YNMjNGQGhjuPnet6Ogou8o2IMTHQSQap2kaACDN/65eAAD4hIABEEnAAIgkYABE\nEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAA\niCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABE+ufqBdjTOI7zL6ZpunZJADhatQls\nmqZpmpaSAVBVqYAZvAD6UeoS4rC6inj3LwF4UfsjQbWAzVt83a32X4MnxnG0/BeKXv7ohR8s/9Ui\nzv7rXEKM2NwA7KXOBLa+dyP6xAeAV9QJ2KBbAD3Jvkr7p/TL0ACXiDh41vkMDICuCBgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARCr1jcwAqcbx\nP3/Z/JdJtsAEBnC1uV7T9PO/4aZn3CNgAA1Yj1zGr9cIGEAbTF1vEjCABoyjwetdAgZwqWXwmn8x\njmL2InchAlxn0yrpeocJDOAit7lSr3eYwABOt9w3zxcEDOBcrhPuRMAAzmLw2pWAAZzC4LU3N3EA\nHE+9DmACAziSy4aHETCAwxi8jiRgAAcweB1PwAD2ZvA6hYC9Zv2UaPslsFgODuvv8XKUOIWAvWC9\nR3rOJrBYHw2k63R5ARt/z3emmx1lXM1Jt7/7leXfNk2+swcYBtcJr5cXsOE3TuM43lZq524N9lHg\nL8vs5ez2XHk/yPw8UeM4jvvuQ+ud0t4J3OU09wqRE9jwYPwa7g1nm57tP6IBfVrObpcPyJMPLzuf\n+p8iL2DzVn5Sr1f+5ns2Q1jyPgrsr8TtG5tDZUTP8gI2PGjSo5lsr//k+r+UvqcCu3E0uE5YwOaT\ngs2NiHO6pml6coMiwM6cy14tLGDPLxKe1K35iqIdF+BSeXchAlzPWWwDBOwjfuADeqZebRAwACIJ\n2KcMYdAn41czBOwLGgZwHQEDeJnxqyUC9h1DGMBFBAzgNcavxgjY1wxh0AP1ao+A7UHDAE4nYAB/\nMX41ScB2YggDOJeAATxl/GqVgO3HEAZwIgEDeMz41TAB25UhDCpRr7YJ2N40DOAUAgZwj/GreQJ2\nAEMYwPEEDOCG8SuBgB3DEAa51CuEgAEQScAOYwiDRMavHAJ2JA0DOIyAAfwyfkURsIMZwgCOIWAA\nwzAYv/II2PEMYdA+9QokYKfQMIC9CRjQPeNXJgE7iyEMYFcCBvTN+BVLwE5kCIPWqFcyAbuCjMHl\nvA3z/XP1AvRnPuNb3jzO/uBk87vPFZF8JrDTzW+baZIuuMDm3SdjyUpNYOPvjji12QZX26Ed3o/5\nqk1g0zRN0zS2eUq1nOs56QP4WqmANTp4Aa1Zxi+nkslKXUKcjeO4LtlmGrs+csvyeOfA+ZbrH26k\n+q9GL1w9VSpg8wuwSdT1xVpbL4xL8HAVb70bm0NlRM9KXUIcWsvVcz4JA/hCnQlsPl9o/UZE4Fqu\nfBRSJ2CRxVp+JgyAN1W7hAjwkPPFWgTsaj4JA/iIgAF9MH6VI2ANMIQBvE/AgA4YvyoSsDYYwgDe\nJGBAdcavogSsGYYwgHcIGFCa8asuAWuJIQzgZQIG1GX8Kk3AGmMIA3iNgAFFGb+qE7D2GMIAXiBg\nQEXGrw4IWJMMYQB/ETCgHONXHwSsVYYwgKcEDKjF+NUNAWuYIQzepV49ETAAIglY2wxh8DrjV2cE\nDIBIAtY8Qxi8wvjVHwEDIJKAJTCEwXPGry4JGACRBCyEIQweMX71SsAAiCRgOQxhcMv41TEBAyCS\ngEUxhMGa8atvAgZAJAFLYwiDmfGrewIGQCQBC2QIA+MXw/DP1QvwoXEcp5vdd1wd1m9/t5S5YbXX\nEeCpvICNT4eP4t0CBuMXP/IuIU7T9KRS4zg+L1wdLiQCfcubwJ6b27a+wLjpWbURbVm7YusFt+zt\nR0o89S8VsLtxqlasxXwVZbmW4qIKta139aJ7+7WB3hwqI3qWdwnxkYjNDXxic6JW8fr5smYVV+4o\nFQI2p2uapvFX2anrlp0d8m1GSm/rF6VeQlwnavl1R92CDhW9csjHKkxg/ZpP0pytUd5mJy9aMu/j\nd6VOYPy8pd2XRVeWk7Zy1lEuGuj9CViy9T5ul6ew0rv3+g6VWd113ZmAVeHhUpBmM096+75LwIC2\nFT0zK7pap3ITRyHu5oAEdX8U+2wmMKBh5Y705VboSiawWgxhVFLuYF9uhS4mYOVoGLTHZcMjuIQI\nNKnQ8b7QqrTFBFaRIQzaYPA6lAkMaE+Jo36JlWiaCawoQxi58g/8Bq9zCFhdGgZXWL7Zi6O5hPiS\ndQjsl3CgqMllc4rogW4nE7C/bfbIpB3U+4ksUbvr5kmGLhuezyXE6lxIhMP4ApRrCdgfNl/SAxwl\nNgKxCx5PwP5QYYCpsA7QKPW6kIC9ZLnYLQRwiMAOrA8IPgC7hJs4/lbhy1LdzUHLkndO760LmcBe\nNf9gR/AQFrzo0KKlW+p1FQF7mxDAnswvfErAeqK9sBPZbYGAfUIIYB86wBcErDPaSzti6xW74NUI\n2IeCQxC86AD/EjDgCrFTTOyCFyRgnwueZIIXHa6kXk0RsK8IAXxCB9iDgPVKe7lKbL1iF7wsAftW\ncAiCFx1AwIAzxU4xsQtemYDtIHiSCV50AokAuxIwgD8ob5sEbB/Bk0zwohMlNgKxC15fwe8DG8dx\numJ3C/5aoE3DIteBVjk94jClAjZ6q3xm+cLp5S81jF0s+1Lse9O7oWWlLiFO03TJ7LVagNj3afCi\n06p1vbK/DZZGlZrA7tqMZdcWDgjS1fiVeAWrfsBOLlb8J2GRi07DYneq2AX/0OZQGdGzUpcQ+ZaL\nPOxoc0rUWxA4Xv0J7Hx5k4xucYT5bbDsWklvibS3cK8KBqyFT7kiG7YWtvS0Z9mF7EgcxiVE7jGT\n8Y3wE6Dwxe+IgB0lPgHxK8BFHP45i4DxmIbRH/0NImAHqnD8r7AOnCj88B+++N0RsGNVOP5XWAdO\n4fDPuQQM2EN+vfLXoDsCdrgKA0yFdeBIjv1cQcB4jYbxSIl6lViJ7gjYGYoc/IusBrty4Oc6AsY7\nNIy1KvWqsh7dEbCT1Dny11kTvlPlqF9lPXokYOepc+SvsyZAMAHjIxrWuSpjS5X16FTBp9G3bHPY\nz37n5D1yn53kv+7OvmowgZ1qfttM08//4t9FFdaBN5Wo1/wGJJ2Anc0xn2BV6jWsMub9mEvAzpP/\n3r/HAaAfNfdgggnYeZZDfbVjfrX14Z5a9aq1Nv0SsMuUOuZrWG2FjvebXbXQmvXIXYinKnUX4oab\nEqsq+rIul0PIJWBnW94wBSeW2znM4SFU2fOsHxXXqUcCdpnKE8u60jXXsLT1q1buPMsuWYnPwNib\nz8Oiberl1aRhAnalsgeHsjdcdqPoB0TGr2IE7GJlj/DzijlgJPKqEULA2NsS5OWBB6S4fbEKxazQ\nqvDDTRzXK3U3xzJ4LX85OHKEuHvvhheOhgkYe7s95K0/EqNBt69OuVfKSVRJLiE2oewnYQuPTW2W\nZ7MTywTGiYxiTenmtTB+VWUCa0Uv84lRrBEGL/KZwLiCUexCnW1541dhJrCG9DWZGMUuYfCiEBMY\nlzKKnabL7Wz8qs0E1pYeZxKj2AkMXlRkAqMNRrGDdLxVjV/l5QVs/D1Vn272zXF1Fn/7uylKPZjj\nLR7bsTsbk9LCAjaO41Km9a8Xud3ix5PLiV7cJ+5utI63mHb3ICxgf5qHsPSM9TuEze6OYj4ke+Lu\ntup3B6IX1QI2p2szqN3+ATI4EL/L17ANw2D8+sgYuM+UCtjdOIUWq/chbOYLWV4n9nxnc6iM6Fmd\ngN39SIwKloz5jo9bd7dJwqHnOM78+hEWsGmabu9CnNN197eiGcK2B2LfVrW43QKbi4fdbhl6Ehaw\n4V6clr9To1v8ePJqdluyJ+vrhxCGYbABOpMXsK4Ywv7QScleX7uSqw8PCBgl/Fmy9dXIpo7yT+JU\nu8oHcLbXGwFrnSHsPbefCd39rUY26HpJll/rFrxGwCiq/RjcdtT9F19o57SE0whYAEPY525vwZ99\ns0Ef3Rv5wb9kc2HTywzvEDD6sPkxqT+fVfEoJK8/4OrPf79c7ce27JOAZXC4+9bt5nu+NZ+XaX19\n8ptHD9/9DAx4jYDRgbn/b30e9nwCW39Y9XF42v+ULoT0d0vAYhjCvrLvhttxbPKKwqf+d/UCQJT1\nMNf3Iwcb4ayuZyawJIawJnhoE7TBBAYfUa8GOIvonICF6f6rCgF+CBgQyfiFgOUxhAEMAgYkMn4x\nCFgoQxiAgAFhjF/MBCyVIYwO2edZ84PMqe5+HQdUtTx+UsNYmMCCzcWaJumiuPma4bKfyxgzAQNi\n+PSLNQHL5lQU6JaAARk23yQKbuIItryNvZ+pbbnS4Os/WROwVJs3sIZRnmix4RIi0Dr3bnCXgBXh\nbg6gNwJWh4ZRkvGLRwQMaJd68YSAlWIIA/ohYECjjF88J2DVGMKoQb34k4AVpGFADwQMaI7xi1cI\nWE2GMKC8Uo+SGn+P2ZOTN4hl/OJFdQI2juPSrfWvuzUPYd1vBsLYaXmdS4iVuZAIFFZnAntk/O8h\n3GQGzTJ+XWgMPNutH7DOi+VCIvCKzaEyomcuIQJNcKbFu+pMYNM0uQvxLkMY7bOL8oE6ARt06zEN\nA+pxCRG4mLMrPiNgvXBLPW1SLz4mYABEErCOGMJojfGLbwhYXzQMKEPAgGsYv/iSgHXHEEYL1Ivv\nCViPNAwoQMCAsxm/2IWAdcoQBqQTMOBUxi/2ImD9MoRxPvViRwLWNQ0DcgkYcBLjF/sSsN4ZwoBQ\nAgYcaz5DMn6xu1JfaMnHNkOYAw27mPcrUz4HMYExDL/FmibpYjfzyLXeo2SMfQkYcKzbksEuBIwf\nrvNwBB99cRwB418axr7W9bJrsTs3cTAMDi4cYD4fWu9aRjH2JWDcOay47MP37EUczSVE7nAtkS+p\nFycQMO7TMD6mXpxDwIA9qRenETAeMoTxLvXiTALGMxrG69SLkwkYf9AwXqFenE/A+JuG8Zx6cQkB\n4yUaxiPqxVUEjFdpGNAUAeMNGsaG8YsLCRjv0TAW6sW1BAz4hHpxOQHjbYYw1IsW5D2Nfvw9dk43\nb6BxdVi9/V12NDfMNu6Tl55GhAVsHMelTOtfL3TrNBrWJy867ah2CXEcx9HlrbO4ltgb9aIpYRPY\nn+YJbDOo3f4B9mIO64cXurbEU/+mA/Zue+7+AcU6moZBAZtDZUTPmg7YW+25+5EY59Cw8ry+NKjp\ngN2apun2LsQ5XXd/i9NoWGFeWdoUFrDhXpyWv6Nb11o3zCEvnZeS9uUFjMbNY/Byg6JjXxyvICmq\n3UZPC+ZDngNfonneWr92CZ/l0ykBY2d+OKyG25JBa1xCZH8als7nXkQwgXEIJ++51vVyIkLLTGDs\nbHPIcyNAkDld4/ifF9FrR7MEjD09Oti5JNW49XmGV4oUAsYZ3JPdMqcXhBIwTuKnYhvkrIJoAsap\njGLtcDJBOgHjbEaxyzmHoAYB4xpGsas4daAMAeMyRrGTOWOgGAHjYkaxczhRoB4B43pGsUM5P6Aq\nAaMVRrEjOC2gMAGjIZtRzAON3rJsrvXWs90oTMBozvph9uunyjoWP3H7BF6bi/I8jZ4WOfi+xfPj\n6ZOA0a71k9F9x9if5m3li2zoh0uINO32wzA25o3jEisdEjACzOPX3c/GurXeFLffwWb70AMBo1Gb\n4/Lmq6q6LdndFd/Mqb1tE7olYLTryYG4t5K9spq1twDcEjCy1S5ZyZWCvQgYRTwv2eZTokZ68Gip\ndAteIWBUc/cnop7E7CqbWy2WHxgYdAteI2CUtS7Z+uEUyz2Nn9n9AVeemAWfETC6cPvzZLdz2Ovx\n+OABV3fHvs0znxoZDSGFgNGR59/b8rwff37nyyv/+ObPN3hhE4IIGPW9+Biq57PU8m+4/anhV/7x\nV5bK9UN4i4DRhe/bcHcC+/KZF4oF3/AwX3iPa33QCBMYvGH9gPzBCAWXEjB4j2hBI1xCbNoYfrnK\n8l8oeuEHy88LCgbMfgPQg1KXEKULoB+lJrBpmiYfUAD0Yax3xB/Hf1fKTAbwmfbrkHoJcVOmRxu6\n/RcAgM+kBkyZADpX6jMwAPpR8DMwAHpgAgMgUupnYK9YbvSInjLXN1Vmid7+0Qu/KLDzDJkvQfT+\n8+Itci0oG7DNzfQtvwaPFPgZgHmzh27/6IUf8vef0M0+/G753P1nvcCN70UuIbYr/eeyLfyFEo+b\nG+M4Nn70fG5e/uhXof3lLzuB0Yj23wOPRB89C8idYIbVCVDo8qcwgXGU9DPQ3Al4Tu/6/+OEbvlK\nIt68AsaB2n8D3BV60F9Mv4bMlyB9+3OagMZ+LPpGoEXEedCtoBuZ7rLzXCt9+xdY/oglz1hKANhw\nCRGASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRg\nAEQSMAAiCRgAkQQMgEgCBkCk/wNa4rq1w/jfzgAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we shift so that one value of $u$ agrees with $v$, then we hope everything will line up:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "u2 = u + v(1)-u(1);\n",
      "norm(u2 - v, inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "text": [
        "ans =\n",
        "     1.567908380906147e-08\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Python versions:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We repeat these examples in Python. The codes are essentially identical, with some changes from Matlab to Python notation.\n",
      "\n",
      "First illustrate how to compute the second derivative of periodic function.\n",
      "Start with $$u = \\exp(\\cos(x)),$$ and check that the numerical approximation agrees well with $$u''(x) = (\\sin^2(x) - \\cos(x))  \\exp(\\cos(x))$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import fft,ifft\n",
      "N = 16;\n",
      "x = linspace(2*pi/N,2*pi,N)\n",
      "\n",
      "ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0)));   # i * wave number vector (matlab ordering)\n",
      "ik2 = ik*ik;          # multiplication factor for second derivative\n",
      "\n",
      "u = exp(cos(x))\n",
      "u_hat = fft(u)\n",
      "v_hat = ik2 * u_hat\n",
      "v = real(ifft(v_hat))     # imaginary part should be at machine precision level\n",
      "\n",
      "error = v - (sin(x)**2 - cos(x)) * exp(cos(x))\n",
      "norm(error,inf)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "3.9095215287332508e-07"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's solve the boundary value problem\n",
      "$$u''(x) = f(x)$$\n",
      "on $0 \\leq x \\leq 2\\pi$ with periodic boundary conditions and the constraint $\\int_0^{2\\pi} u(x) dx = 0$.\n",
      "\n",
      "Use $f(x) = (\\sin^2(x) - \\cos(x))  \\exp(\\cos(x))$ so the solution should be $u(x) = \\exp(\\cos(x)) + C$, where the constant is chosen so the integral constraint is satisfied."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "N = 16;\n",
      "x = linspace(2*pi/N,2*pi,N)\n",
      "f = (sin(x)**2 - cos(x)) * exp(cos(x))\n",
      "f_hat = fft(f)\n",
      "ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0)));   # i * wave number vector (matlab ordering)\n",
      "ik2 = ik*ik;          # multiplication factor for second derivative\n",
      "ik2inverse = where(ik2 != 0, 1./ik2, 0.)\n",
      "u_hat = ik2inverse * f_hat;\n",
      "u = real(ifft(u_hat))\n",
      "\n",
      "plot(x,u,'b-o')\n",
      "v = exp(cos(x));\n",
      "plot(x,v,'r-o')\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:8: RuntimeWarning: divide by zero encountered in divide\n",
        "-c:8: RuntimeWarning: invalid value encountered in divide\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10841f390>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD9CAYAAACyYrxEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xlc1NX+P/AXAoGYWy6gYGJjubC7YfYTR5NBRVGv1VVb\nNJevXYtRu7dSkCteRbtpqaCleXPrVrdb1wUZF6gcxiWXEsFMDUdKcMHUTFBZZji/P04iI4POfuYz\n834+HvNQhg/MC/jMez5zVg/GGAMhhBCX1Uh0AEIIIfZFhZ4QQlwcFXpCCHFxVOgJIcTFUaEnhBAX\nR4WeEEJcnFWFvqKiAtHR0YiMjETfvn2xbNkyo8fNmTMH4eHh6Nu3L06dOmXNQxJCCDGTlzVf7Ovr\niz179sDPzw+VlZXo2bMnRowYgc6dO9ces2PHDuTn56OgoACHDh3CxIkTcfDgQauDE0IIMY3VTTd+\nfn4AgPLycuh0Ovj4+Bh8PjMzExMmTAAAREdH4/r16ygtLbX2YQkhhJjIqit6AKipqUFUVBROnDiB\n5cuXo0OHDgafP3/+vMF9QUFBKCkpgb+/f+19Hh4e1sYghBC3ZMriBlZf0Tdq1Aj5+fk4c+YM3n//\nfeTl5T0wiLHCzhiT7G3evHnCM7hjdsov/kb5xd5MrtPml3bjgoODMWzYMOTm5hrcHxgYiOLi4tqP\nS0pKEBgYaKuHJYQQ8gBWFforV67g+vXrAICrV69i586dCAsLMzgmISEBmzZtAgAcPHgQLVq0MGi2\nIYQQYl9WtdFfvHgREyZMgF6vR0BAAF5//XU8/fTTWLNmDQBg2rRpGDZsGDQaDcLCwtCkSROsX7/e\nJsGdiVwuFx3BYlLODlB+0Si/NHgwcxp67BXCw8Os9iZCCCGm106aGUsIIS6OCj0hhLg4KvSEEOLi\nqNATQoiLo0JPCCEujgo9IYS4OKvXuiGEEHegUamQnZ4Or8pK6Hx8oFAqERMfLzqWSajQE0LIA2hU\nKuyeMQNpWm3tfcl//F8KxZ6abggh5AGy09MNijwApGm1yMnIEJTIPFToCSHkAbwqK43e71lR4eAk\nlqFCTwghD6C7Z0OlO/S+vg5OYhkq9IQQ8gCK8eOR/PDDBvclyWSITUwUlMg8tKgZIYQ0pKoKWLYM\nWLIEmqefRs61a/Csrobe1xexiYnCO2JNrZ006oYQQozZtQuYMQN4/HHg4EHEdO6MGNGZLESFnhBC\n6ioqAmbNAn74AVi+HBg+XHQiq1EbPSGEAMCtW8C8eUDv3kCfPrzQu0CRB+iKnhDi7hgDtmwBXn8d\niI4G8vKADh1Ep7Iplyj0Up6aTAixnwfWhlOnAKUSuHABWLcOGDRIXFg7knyhl/rUZEKIfdy3NvTv\nDyxYAGzYACQnA6++Cnh7C0pqf1a10RcXF2PgwIEICQmBXC7Hhg0b6h2jVqvRvHlzREVFISoqCgsX\nLrTmIeuR+tRkQoh9NFgb5swBunUDrlzh7fAzZ7p0kQesvKL39vbGsmXLEBkZiStXriA0NBTR0dHo\n1q2bwXEDBgxAZmamVUEbIvWpyYQQ+2iwNhQVATk5QN++Dk4kjlVX9AEBAYiMjAQAtG7dGr1798aF\nCxfqHWfPyVBSn5pMCLGPBmtDv35uVeQBG7bRnzlzBidOnEDfe36BHh4eOHDgAEJCQtCxY0csXboU\n3bt3r/f1qamptf+Xy+WQy+UmPa5CqUSyVmvwFi3poYcw5JVXLPo5CCGuwWhtkMkwRKkUmMo6arUa\narXa7K+zyRII5eXlkMvlSElJwciRIw0+V1ZWBk9PT3h7e2Pjxo14++23cebMGcMQVi6BoFGpkJOR\nAc+KCj41+do1xMTGAmlpFn9PQoj0aVQq5EyfzmtDVJRTLFtgS6bWTqsLfXV1NYYPH46hQ4di5syZ\n9z2WMYbWrVujsLAQjzzyiNlhTXbpEhAZCWzd6nZv0QghdXz1FfDyy0BBAdCypeg0Nmdq7bSqjZ4x\nhsmTJyMkJKTBIl9aWlobZPv27WjcuLFBkbeLgABg1SrgpZeAmzft+1iEEOd0/TowaRLw0UcuWeTN\nYdUV/b59+xATE4Pw8HB4eHgAABYtWoRz584BAKZNm4ZVq1bhgw8+gJeXF8LDwzFjxgz07NnTMIS9\nVq984QX+B6ahloS4n5deAh5+GHj/fdFJ7MZhTTe2YLdC/9tvQHg4sH49MHiw7b8/IcQ5bd4MvPUW\ncOwY0KSJ6DR2Q4X+juxsYMoU3kbXooV9HoMQ4jxKS4GICF7s+/UTncauqNDX9eqrQFkZsGmT/R6D\nECIeY8Do0Xzm6+LFotPYnUM6YyXjnXeAb7/lr/CEENe1cSNfT77OvBziLlf0AHDgAPCnPwH5+YC/\nv30fixDieL/8AvTqxYdURkSITuMQ1HRjzJw5wMmTfO3pP0YJEUJcQE0NH3ARG8uf526Cmm6MSU3l\nb+s2bhSdhBBiSytXAhUVwBtviE7ilNzrih7go28GDwaOHAE6dnTMYxJC7OfUKeD//T/eD/f446LT\nOBRd0TckPJxvGfbyy/ztHiFEunQ6PjHqH/9wuyJvDvcr9AB/e1dZyd/uEUKka/FiPj/mL38RncSp\nuV/TzR1nzgBPPgns3Qt07erYxyaEWO/oUWDIEP5vUJDoNEJQ082DdO7M3+5NmMDf/hFCpKOiAnjx\nReC999y2yJvDfQs9ALzyCn/b9/bbopMQQsyRksJnvz7/vOgkkuC+TTd3lJQAPXoAu3bxfwkhzk2j\nAcaO5ZMf27QRnUYoaroxVVAQsGwZ77mnDcUJcW5lZcDEicDq1W5f5M1BV/QAXwjpueeA4GBgyRJx\nOQgh9/d//wfo9XwzEUJLIJjtyhU+xv7zz4H+/cVmIYTUt2MHMH06n/TYrJnoNE6Bmm7M1bo1sGYN\nH4VTViY6DSGkrqtXgalTgQ0bqMhbgK7o76GJjUX2qVPwksmg8/GBQql0qV3jCZECjUqF7PR0eFVW\n8ufh7duI6dmT96eRWqbWTi8HZJEMjUqF3Vot0kpK+GgcAMlaLQBQsSfEQTQqFXbPmIG0P557AJDs\n7Q3MnIkYgbmkzKqmm+LiYgwcOBAhISGQy+XYsGGD0ePmzJmD8PBw9O3bF6dOnbLmIe0qOz0daUVF\nBvelabXIoc3FCXGY7PR0gyIPAGnV1cj58ENBiaTPqit6b29vLFu2DJGRkbhy5QpCQ0MRHR2Nbt26\n1R6zY8cO5Ofno6CgAIcOHcLEiRNx8OBBq4Pbg1dlpdH7PWnYJSEOQ89D27Oq0AcEBCAgIAAA0Lp1\na/Tu3RsXLlwwKPSZmZmYMGECACA6OhrXr19HaWkp/O/Z5Sm1ztZfcrkccrncmmgW0fn4GL1f7+vr\n4CSEuC96HjZMrVZDrVab/4XMRgoLC1mnTp1YeXm5wf3Dhw9n+/fvr/346aefZt99953BMTaMYZXc\nrCyWJJMxxkfWMwawOd7eLHfLFtHRCHEbuVlZLOmxxwyfhzIZy83KEh3N6ZhaO23SGVteXo6xY8di\n2bJlaNKkibEXE4OPPZx0G787Ha4pGRnwrKiA3tcXQ65eRUxhoeBkhLiPmPh4YNMmpNy6Bc8uXfjz\nMDGRBkRYwepCX11djTFjxuCFF17AyJEj630+MDAQxcXFtR+XlJQgMDDQ2oe1m5j4eMMT6swZoG9f\nYNw4WiWPEEcoLETM118jJj8fcOJaISVWjbphjGHy5MkICQnBzJkzjR6TkJCATZs2AQAOHjyIFi1a\n1Gufd2qdO/PZeH/7m+gkhLg+xgClEpg9m4q8DVk1YWrfvn2IiYlBeHh4bXPMokWLcO7cOQDAtGnT\nAACzZ8+GSqVCkyZNsH79eoPOWsC5JkwZdesWEBLC19cYNEh0GkJc19atQFISX5nS21t0GqdHa93Y\n2p0T8Ngx4KGHRKchxPXcugV07w6sW0cXVCaitW5sbeRIvrplerroJIS4psWLeX8YFXmboyt6c9zp\nmKVOIkJsq7CQ7+F87BgNejADNd3Yy9y5gFYLfPaZ6CSEuAbGgPh4YOBA4I03RKeRFGq6sZekJODb\nb4E9e0QnIcQ1ZGYCRUXAjBmik7gsKvTm8vPjS6W+9hpQXS06DSHSdusWMHMmsHIlDXKwIyr0lhg1\nCujQgTpmCbHW228DffoATz8tOolLozZ6S93pPCooANq3F52GEOm5M7iBOmAtRp2xjpCczNsWP/1U\ndBJCpIUxYPhwYMAA4M03RaeRLOqMdYSkJGD/fsCSZUMJcWfbtwNnz/L2eWJ3VOit0aQJ75h99VXq\nmCXEVLdv8xE2GRnUAesgVOitNXo0b1+k7QYJMc3bbwO9ewODB4tO4jaojd4WfvoJ6NePOmYJeRCt\nFoiOBvLy+Mg1YhXqjHW0pCTgl1+ATz4RnYQQ5zV8ONC/P/DWW6KTuAQq9I528ybQrRvw8cd8JAEh\nxND27XyJg4ICapu3ERp142jUMUtIw6gDVigq9Lb0pz/xNvqVK0UnIcS5/POfQM+eQGys6CRuiZpu\nbO30aeCpp4Djx4F27USnIUS8s2f5MgfUAWtz1EYv0pw5QHEx8O9/i05CiHgjRvCLn9mzRSdxOQ5p\no580aRL8/f0RFhZm9PNqtRrNmzdHVFQUoqKisHDhQmseTjrmzgU0Gn4jxJ1lZfHhx6+/LjqJW7Oq\n0L/88svYtWvXfY8ZMGAA8vLykJeXh7lz51rzcNLRpAnw7rvUMUvcG3XAOg0va764f//++Pnnn+97\njEs1yZjjmWegSUtDdlgYvAICoPPxgUKpREx8vOhkhNiFRqVCdno6vCor+fneujVioqIAhUJ0NLdn\nVaF/EA8PDxw4cAAhISHo2LEjli5diu7duxs9NjU1tfb/crkccrncntHsTrNjB3Zfu4a04mLeQQsg\nWasFACr2xOVoVCrsnjEDaX+c4wCQ3KgR8NFHiBGYy5ZUKg3S07NRWekFHx8dlEoF4uMd+9Op1Wqo\nLVlEkVmpqKiIhYaGGv3cjRs32M2bN1lVVRVbu3Ytk8lkRo+zQQynk6xQMMYXYzW4zY2LEx2NEJtz\n9fM9KyuXyWRJBj+eTJbEsrJyheYytXbadRx906ZN4efnB29vb0yePBm//fYbrl27Zs+HdBpelZVG\n7/esqHBwEkLsz9XP9/T0bGi1aQb3abVpyMjIEZTIPHYt9KWlpbVt9Nu3b0fjxo3xyCOP2PMhnYbO\nx8fo/XpfXwcnIcT+XP18r6w03spdUeHp4CSWsarQjxs3Dv369cPp06fRoUMHrFu3DmvWrMGaNWsA\nAF9++SXCwsIQGRmJL7/8Etu2bbNJaClQKJVIlskM7kt6+GHEJiYKSkSI/Rg932UylznffXx0Ru/3\n9dU7OIllaMKUHWlUKuRkZMCzogJ6b2/EnjyJmHfeAcaPFx2NEJvTvPoqcj7+GJ5RUdA3bozYxESX\nGHig1wNKpQarV+9GTc3d5huZLAkrVgxxeIdsXTQz1hnl5fGhZocOAY89JjoNIbZz4gQglwP79gFd\nuohOYzP79wOvvQY0awY8+6wGWVk5qKjwhK+vHomJsUKLPECF3nktWwb897981qy3t+g0hFivooKv\nZTNzJjBpkug0NnHxIl8y/5tvgCVLgLFjAQ8P0anqo2WKndWMGUDz5sD8+aKTEGIbb7wBdO0KvPyy\n6CRWq6oCli4FwsKAwEDg1Clg3DjnLPLmsOuEKWJEo0bAxo1AVBTfM1PiE8OIm9u+nd+OHZN8NczJ\nAZRKIDgYOHAAeOIJ0Ylsh5puRNm1C5g6lT9BWrUSnYYQ8124APToAfzvf3x1Son65Re+5lpeHrB8\nOV9sUyqvWW7VRu8MU5Mt8vrrQFERsHmzdM4sQgCgpoYPLOjfH5g3T3SaBt2vNty+zdvfV6zg3Qt/\n+xvQuLHgwGYytXZKvulGpdJgxozdBrPWtNpkAHD+Yr94MdC3L7BmDfDKK6LTEGK6pUuBykogOVl0\nkgY1VBsYA/T6GMyaxd+QHD0KdOwoMKgj2G7VBctZE0OhSDa2xAaLi5trw4R2dPIkY61aMXbihOgk\nhJjm8GHG2rRh7OefRSe5r4ZqQ+vWc1nXroxlZ4tOaD1Ta6fkr+ilPjUZXbsCb7/Nx28dPgy4yJRx\n4qLKyvgwlFWrnP4yuKHa0Ly5J/Lz3WuJfMkPr5T61GQAwOTJfJLJm2+KTkLI/b32Gh8p9uyzopM8\nUEO1oXNnvVsVecAFCr1SqYBMZthO2KhREp59VkK7zXt4AB9+CGzbxrdeI8QZffopn9W9YoXoJCYx\nVhtksiQkJkqoNtiIy4y6yci4OzW5fftYaLUx2LOHD1uXjH37gGee4eO82rUTnYaQu4qKgOhoYPdu\nPgdEIrZt0+DFF3PQrp0nOnVyjmULbMmthlfeS68HBg4ERo2S4J7Eqal8gY3duyX2KkVcVnU1EBMD\nPPccMGuW6DRmmTcPOHIEUKlccwSzWxd6ADh7ll+A5OYCDexe6Jx0Ot4GOnIkn1pOiGhz5wLffQfs\n2CGpi48jR4Dhw/kb5PbtRaexD7cv9ACwdi2wejVw8KDE1g/75Regd2/+xOrVS3Qa4s7Uar6sdl4e\n4O8vOo3Jbt/mLUypqXxAm6uiRc0ATJnCz820tAcf61Q6dgRWruTD2MrKRKch7uraNeCll4B16yRV\n5AFgzhwgMtK1i7w5XPqKHuDLjUZG8sEsvXvb5SHsZ/Jk3uGwYYPoJMTdMAaMGcNX+HrvPdFpzLJn\nD/Dii0B+vusvI0VX9H9o1w5IT+cXJrdvi05jphUrgG+/BT77THQS4m4+/JCPtFm8WHQSs/z+O18t\n+cMPXb/Im8Plr+jvGDcOCAjg+35IytGjQFwcnzXbqZPoNMQd/PgjMGAAsHcvn7ktIZMm8f64P7at\ndnkOWdRs0qRJUKlUaNu2LY4fP270mDlz5kClUsHPzw8bNmxAV0EnzqpVQHg4kJDAh15KRo8ewOzZ\n0AwdiuxHH4VXVRV0Pj5QKJUusR8nEUujUiE7PR1elZX8vHrlFcTMm8ev5CVW5DMzed9xfr7oJE7I\nmgV1NBoNO3r0KAsNDTX6eZVKxYYOHcoYY+zgwYMsOjra6HFWxjDZjh2MdezI2O+/O+ThbCY3M5Ml\nNW5ssDJTkkzGcrOyREcjEpablcWSZDLD86p5c5b71FOM1dSIjmeWy5cZa9eOMY1GdBLHMrV2WtVG\n379/f7Rs2bLBz2dmZmLChAkAgOjoaFy/fh2lpaXWPKRVhg7lrSAzZwqLYJHslSuRdk8HQ5pWi5yM\nDEGJiCvITk9HmlZrcF/a778jx8dHUrOLGOOrfD//PF8en9Rn19Urz58/jw4dOtR+HBQUhJKSEvgb\nGaqVmppa+3+5XA65nbbYe/ddICKCv81LSLDLQ9icV2Wl0fs9KyocnIS4kgbPK72EFgQE8MknwOnT\n/F9Xp1aroVarzf46uy9TzO7pKPBo4EqhbqG3p4cf5qMV//xn4MkngTZtHPKwVtH5+Bi9X09LGhMr\nuMJ5VVzMlznZvds9Vvi+9yJ4/vz5Jn2dXYdXBgYGori4uPbjkpISBAYG2vMhTdK/P/DCC/ztnvgx\nRw+mUCqRLJMZ3Jfk7Y1Ymg1CrKBQKpH86KMG9yXJZIhNTBSUyDw1NXyUjVIpqXXWhLDrFX1CQgJW\nrlyJsWPH4uDBg2jRooXRZhsR/vEPvrrAJ5/wou/M7oyuScnIgGdFBfS+vhjSvj1i5s/nr1r3vAgQ\nYoqY4GCgvBwp3bvDs00bfl4lJkpmNNcHHwA3bgCzZ4tO4vysGkc/btw45Obm4sqVK/D398f8+fNR\nXV0NAJg2bRoAYPbs2VCpVGjSpAnWr1+Pbt261Q/hgHH0xuTl8c7Zo0eBoCCHP7z11qwBFizg71tD\nQkSnIVLy/fdAfDyf9Tp+vOg0Ziss5E2v+/fzPXvcFS1qZqKFCwGNhtdKCQ00uOvTT3kjZVYWLYBG\nTLN3L1/eYO1avkqqxOh0/I3suHG82cad0RIIJpo9m0+b/uAD0UksNH48v7IfNoy/YhFyP7t28SL/\n6aeSLPIAsGQJ4OfHdzUkpnH7K3qAD8166im+rMzjjwuLYZ2vvuKXOB9/DAwZIjoNcUb/+x8wfTqw\nZQvQr5/oNBbJzwcGD+YtT/f0I7sluqI3Q5cufCeaCRP420JJGjyY7zk7YQLw5Zei0xBns3EjkJjI\n2yglWuQrK/mqlEuXUpE3FxX6P7z6Kn87uGSJ6CRW6NePP5ETE2lpY3LXypVASgrwzTd8zW6JmjcP\neOwxvhItMQ813dRRXAyEhmoQEpKNhx7ygo+PDkqlQnqbCZ86BSgUfCtCiYyJJnayeDHw0Ue8aS84\nWHQak6lUGqSnZ6Oykj8Pn35agWXLYpCfD7RtKzqd83DI6pWupqBAA1/f3fj227tbUmm1yQAgrWLf\ntSsfWTF4MN+has4ciQ4pIhZjjP/ds7L4udCunehEJlOpNJgxYze02rvPw2++ScabbwJt20roeehE\nqOmmjvT0bFy+bLjvoFabhoyMHEGJrNCxIx+F89lnfGiRE7xjIg5SU8PbIr/+GsjNlVSRB/jzsG6R\nBwCdLg3ffy/B56GToEJfR2Wl8Tc4FRWeDk5iI+3a8QW69+zhoy1qakQnIvam0wETJwInTvBCL8Ft\nllzueegEqNDX4eNjfMiNr6+0VvMz0KoVb5/98UeJDysiD1RZCTz3HPDrr8DOnUCzZqITWcQln4eC\nUaGvQ6lUQCZLNrjvkUeSkJgYKyiRjTRrxp/4V68Czz4LzdatmBsXh1S5HHPj4qBRqUQnJGbQqFT1\n/343bwIjRgCennyYrZ+f6JgWM/Y8lMlc4HkoEHXG1nGnwzUjIwUVFZ5gTI+jR4cgLMwFOoD8/ICt\nW6EZNAi7x4832Mgk+Y/NJ6SymJU706hU2D1jhsGGIcmFhYCvL2Kio/myBl7SflrHx8eguBhQKlPQ\nq5cnmjXTIzFxiLQGRDgZGl75APPnA8ePu84cpLkKBRbm1O/USomLw4JduwQkIuaYGxeHhdnZ9e5P\nefRRLCgqAhpJ/006Y3x0cHy89HaDczSaGWsjb77JV7k08tySJK+qKqP3025V0tDgrlCdOrlEkQf4\nSg2XLtFaNrbkGmeGHTVuDKxYwecdNfAck5QGdxVykSLh6lxhV6j7KS/ni7GuWiX5FiinQs9uEwwf\nztfDWbZMdBLrGd2tqnlzxB45wtd/oCt753XrFhRt2iD5nhdlKe0K9SBpaUBMDL8R26E2ehOdPQv0\n6cObcersdy5JGpUKOXV2q4pNTETM44/zdqr8fOCf/wSefZZm0zqLmhrg3/8GkpOBfv2gUSiQ88UX\nhn8/F+hIv7OK7PHjkpvjJQxtPGIHqal8HsoXX4hOYkd79vD3zo0b892H+vYVnci95ebyv4e3N/97\nSHTlyQdhjO/2NnQoMGuW6DTSQZ2xdvDWW3wdbCODVlzHwIHAd98B//d/wDPP8DXuf/5ZdCr3U1gI\njB7NJ7m98QbfLMFFizwAbN4MXLhAHbD2QoXeDHU7ZhsYvOIaPD35NPrTp/kCaT178gWybtwQncz1\nXbvGL2mffJK/mzp1Chg71qWb0W7e5D/yqlX8jQuxPasLvUajQY8ePRAeHo6MjIx6n1er1WjevDmi\noqIQFRWFhQsXWvuQQo0YAXTu7Bodsw/UpAlfBLyggI9369KFb1v4xzIKRmdokgbd9/dVVQUsX85f\nWCsr+ZIVb70FuMhomvtJS+N7wA4YIDqJC2NW0Ol0TCaTsaKiIlZVVcUiIiLYjz/+aHDMnj172IgR\nI+77fayM4XBnzjDWqhVj586JTuJgR48yNnAgY927s9zUVJYkkzHGm1cZA1iSTMZys7JEp3RKuVlZ\nxn9f27cztmULY507MzZ0KGM//CA6qkOdPs2fS+fPi04iTabWTqsq7IEDB1hcXFztx4sXL2aLFy82\nOGbPnj1s+PDh9w8hsULPGGN//ztjzz4rOoUANTWMbdvGkv38DIrWndvcOucDuStZoTD++2rZkrHQ\nUMZ27xYd0eFqahhTKBhbulR0EukytXZaNSXh/Pnz6FBnrGFQUBAOHTpkcIyHhwcOHDiAkJAQdOzY\nEUuXLkX37t3rfa/U1NTa/8vlcsjlcmui2d3s2UD37nxhyMGDRadxIA8PICEBXr168fXu70EzbI1r\ncEZrmzbAsWO8X8TNbNkClJQASqXoJNKhVquhVqvN/jqrCr2HCR1EPXr0QHFxMby9vbFx40YkJCTg\nzJkz9Y6rW+iloG7HbH4+8NBDohM5lq6BtmN9aSnvxH3iCZfuQDQLY9BVVxv9lL5TJ7cs8nc6YDdu\npA5Yc9x7ETx//nyTvs6qztjAwEAUFxfXflxcXIygoCCDY5o2bQo/Pz94e3tj8uTJ+O2333Dt2jVr\nHtZpjBjBNytevlx0EsczOsPW3x+xjz4KPP00IJPxXY6ysviz2t2UlQFbtwLTpgEdO0JRWIjke9aH\nd6UZreZatIhPjnLyN+4uw6oJUzqdDl26dMHXX3+N9u3bo0+fPvjss8/QrVu32mNKS0vRtm1beHh4\nIDMzE9OnT0dJSYlhCIlMmDLmzBk+Cu7YMeCe1ziXZ3SGbXw8b33+4Qe+Bv7OnXxc/pNP8tkww4YZ\nvdrXqFTITk+HV2UldD4+UCiVTjHb0+RcjPHZdHd+5iNH+IkxdCi/de0KzY4dxn9fbuann/iUgIIC\noH170WmkzeTaaW1ngFqtZpGRkSw0NJStWLGCMcbY6tWr2erVqxljjK1cuZKFhISwiIgI9uKLL7Lv\nvvvO4g4FZ5WSwthzz4lO4cR+/52xzZsZmzqVsaAgxjp1Ymz6dMa2b2esvLzhESkWjuDJzcpiyQoF\nmzdgAEtWKKz6PvfNdePG3Z+rQwfGgoMZ+8tfGMvMZKyszKLHdHU1NYzFxTG2ZInoJK7B1NrpFBVW\n6oX+5k3+HP/qK9FJJKCmhrHjxxl75x0+VPPhh1nyI4/YbASPLV80Ghwp88QTtdlZbCxj773H2MmT\n/Gcj97V5M2PduzNWVSU6iWswtXbSQqA24OfHJ1AlJvImHHfrmDWLhwcQGspvb7wBlJXB68kn+YzQ\ne3gWFgLij6UOAAAUt0lEQVSbNgFt2969tWkDNLBULwBkp6cb7L4EAGlaLVIyMh7cTFJRwfdbvXwZ\nKC2F1y+/GD3M885augMH8kllxCS3bvEO2PXrqQPW0ajQ28jIkcCHH/KROG+8ITqNhDRtCl1gIG/f\nvoder+cLC12+fPf266/8lbVu8W/bFvD3B9q2hdc9/T93eF66BPz3v0BpqeH3q3u7fdvg++kaWPJB\nHxbG164mZlm0iHdbDBwoOon7oUJvIx4eQHo6P5HHjwcCA0Unkg6FUolkrdbgSjxJJsOQFSv4fnJ1\nMQZcv25YoO8U7xMnoLtyxehj6M+d48uO/vGCgKio+i8WzZsbdBIrVCok37M/a5JMhiFuOlLGGoWF\nwOrVfCgycTxaptjGUlL4SJzPPhOdRFoaHMFjwfe5d/PsOy8aln4/GiljHcb46/XAgfRu19ZoPXpB\nbt3iM2bXr6e3qKJQcXYuW7cCSUnUf2UPVOgF2rKFbwaUn0+dTsS93bnwWbcOGDRIdBrXQxuPCDRq\nFNCxI++YJcSdLV4MREdTkReNrujtpLCQTwal2X/EXbnzrHFHMbV20qgbO3n8cWDQIA169sxGly5e\n8PHRQalUID6etrcnrkml0iA9PRuVlfx8v3pVgTffjKEi7wSo0NuJSqXB99/vxqVLabh0id+n1SYD\nABV74nJUKg1mzNgNrTat9j5v72TMnQsAdL6LRm30dpKeno2zZ9MM7tNq05CR4co7ixN3lZ6ebVDk\nAaC6Og2rV9P57gyo0NtJZaXxN0sVFe639jhxfXS+Ozcq9Hbi46Mzer+vr97BSQixPzrfnRsVejtR\nKhWQyZIN7vP2TsLUqbGCEhFiP0qlAp06GZ7vMlkSEhPpfHcGNLzSjlQqDTIyclBR4QlfXz2qqmIh\nk8Vg7VrRyQixLcaA/v01KC7OQadO/HxPTIylgQd2RjNjnVBZGV9L6+23gWeeEZ2GENv56CM+QfDw\nYaCB7YSJHVChd1JHjvAFnr77Dnj0UdFpCLHeqVNA//5Abi5f7oA4Di2B4KR69wb++lfg+ecBnfH+\nK0Iko7KSL8u9YAEVeWdGV/QC1NQACgW/Cpo3T3QaQiz3178CZ88CmzfX2++dOIDDrug1Gg169OiB\n8PBwZGRkGD1mzpw5CA8PR9++fXHq1ClrH1LyGjXiO+R98AGwf7/oNIRYZtcuvmnXv/5FRd7ZWVXo\n9Xo9Jk2ahM2bN+P777/HRx99hJMnTxocs2PHDuTn56OgoAArVqzAxIkTrXlIl9G+PbB2LW/CuX5d\ndBpCzFNaCkyaxC9YWrUSnYY8iFWF/vDhw+jcuTOCg4Ph7e2NsWPHYtu2bQbHZGZmYsKECQCA6Oho\nXL9+HaWlpdY8rMsYMYJvPTptGh+eRogU1NQAEycCL79Mm+tIhVWLmp0/fx4dOnSo/TgoKAiHDh16\n4DElJSXw9/c3OC41NbX2/3K5HHK53JpokrFkCdCnD9+RatIk0WkIebD0dOC334A6T1niIGq1Gmq1\n2uyvs6rQe5jYMHdvZ4Gxr0t107OmcWO+v6xcDjz1FNCli+hEhDQsLw9ISwMOHaLd00S49yJ4/vz5\nJn2dVU03gYGBKC4urv24uLgYQfcsPn3vMSUlJQgMDLTmYV1OaCjwj38A48bx4WqEOKObN/k5unw5\n8NhjotMQc1hV6Hv16oXCwkL8/PPPqKqqwueff46EhASDYxISErBp0yYAwMGDB9GiRYt6zTYE+Mtf\n+ASq5OQHH0uICDNn8mbG558XnYSYy6qmGy8vL6xbtw6jR4+GTqfD1KlT0a1bN6xZswYAMG3aNAwb\nNgwajQZhYWFo0qQJ1q9fb5PgrsbDgw9Ti4oCYmOBuDjRiQi568svgT17eNMNkR6aMOVk9uzhV0zH\njgFt24pOQwhw7hyf0Z2Vxf8lzoPWupGwpCRe6LOy+OQqQkTR6/kQyvh44K23RKch96K1biRs/nzg\n6lWggYnGhDhMWhofXfPGG6KTEGvQFb2T0mqBvn2BnBwgMlJ0GuKO9u8HxowBjh7lM7mJ86EreomT\nyYBly/hwtps3Rach7ub6dd5XtHYtFXlXQFf0Tu7FF/mkqg8/FJ2EuAvG+AVG69bAypWi05D7oSt6\nF7FqFfD113x4GyGOsGEDcOIEX56DuAa6opeAQ4cAhUKDyMhseHh4wcdHB6VSQftxEqupVBqkp2ej\nspKfV2PGKJCcHIM9e/iMbeLcTK2dVk2YIo5x5YoGXl67odGk1d6n1fIptFTsiaVUKg1mzNgNrfbu\neaVWJ2PyZCA0lM4rV0JNNxKQnp6Na9fSDO7TatOQkZEjKBFxBenp2QZFHgCqqtJw9iydV66GCr0E\nVFYaf+NVUeHp4CTEldB55T6o0EuAj4/xXcR9ffUOTkJcCZ1X7oMKvQQolQrIZIbLWjZqlITw8FhB\niYgrUCoVaNPG8LySyZKQmEjnlauhzlgJuNPhmpGRgooKT/j66jFixBAsWhSDxx8Hpk4VHJBI0k8/\n8fOqf/8UNGrEz6vExCHUwe+CaHilhBUW8iWNlUrg9ddFpyFSwRjf6OaTT4CvvuL7IBBpouGVbuDx\nx4G9e4HBg4EbN4B58/i69oQ0hDHgb3/jBX7vXoD2AHIPVOglrkMHQKPhG5XcuAG8+y4Ve2KcXg+8\n8gpw/DigVgMtW4pORByFmm5cxG+/AcOG8dmMq1cDnjRCjtRRXc3XTbp8Gdi2DWjaVHQiYgu01o2b\nadmSL2l89ixfdbCqSnQi4ixu3wb+9Ce+CuqOHVTk3REVehfy8MOASgXcusWf2Ldvi05ERCsr47tD\nPfwwsHkz4OsrOhERweJCX1ZWhlGjRiE8PByjR49GeXm50eOCg4MRHh6OqKgo9OnTx+KgxDS+vsD/\n/gc0a8abcsrKRCcioly7xkdlyWTAv//Nd4oi7sniQr9gwQL069cPBQUF6Nu3LxYuXGj0OA8PD6jV\nauTl5eHw4cMWByWm8/YGPv6Yj8oZPJg/4Yl7KS3le70+9RTfy4D6bNybxYU+MzMTEyZMAABMmDAB\nW7dubfBY6mh1PE9PYM0aoH9/QC4HLl0SnYg4yrlz/O8+ZgywdCmNwiJWDK8sLS2F/x+DcP39/VFa\nWmr0OA8PDwwaNAiNGjXC9OnTMbWBaZypqam1/5fL5ZDL5ZZGI3/w8OCbRzRvDsTE0OQYd3BnEt2M\nGcCsWaLTEFtTq9VQq9Vmf919h1fGxsbikpFLwbS0NEyYMAG//fZb7X2PPPIIrhlpI7h48SLatWuH\nkydPYtiwYdi0aRP69+9vGIKGV9rd8uV8D9qcHOCJJ0SnIfZQUAAMHQrMnw9MmSI6DXEEm8yMzclp\neF1qf39/XLp0CQEBAbh48SLatm1r9Lh27doBALp164bRo0fj8OHD9Qo9sb+ZM/mwOrkcSErSYPv2\nu7sK0W5V0nLvrlBKpQKtW8cgIQFITwf+/GfRCYmzsbjpJiEhARs3bsRbb72FjRs3YtSoUfWOuXXr\nFvR6PZo2bYpff/0VO3bsQHp6ulWBieUmTwZ++onvKlRTQ7tVSZGxXaF++CEZ5eXAp5/GID5eYDji\nvJiFbty4wUaOHMnCwsLYqFGjWFlZGWOMsfPnz7Nhw4YxxhjTarUsIiKCRUREsEGDBrHVq1cb/V5W\nxCBmUiiSGV/xxPAWFzdXdDRigob+fr1709/PHZlaOy2+om/atKnRkTbt27eHSqUCADz22GM4duyY\npQ9B7IB2FZK2hv5+fn709yMNo5mxbqahXYWuXtWjpsbBYYjZHnqIdoUi5qNC72aM7VYVGJgEnS4W\nvXoBe/YICkYeaNcu4NQpBRo3pl2hiHlo9Uo3pFJpkJGRU7tbVWJiLIYNi8EXXwBvvQVERADvvEPD\nMJ3FDz/wNeSLivi8iEaNNFi50vDvRx3p7snU2kmFnhioqOBD9N55h6+C+fe/A61aiU7lnkpL+e9/\n61Zg7ly+ljytV0PqomWKiUV8fYE33wROngR0OqBrVz7RipY9dpzbt4HFi4GQEL7q5KlTQGIiFXli\nOSr0xKg2bYBVq/juVV9/DXTvzpe5pTde9lNTA3z6KX9x/f574NAhvmMY7QRFrEVNN8QkOTnAX/8K\ntGgBvPce0KuX8Rma1FbcsPv9vvbt4xu8M8Z/vzR5nJiCNgcnNhUbC+TlAevXAwkJQJcuGhQV7cYv\nv9AMW1MYm9Gq1Sbj4kVg9+4YHDrEm2vGjQMa0ftsYmN0RU/MVlYGREbOxdmz9fcgiItLwa5dCwSk\ncm5xcXORnV3/9+XllYLU1AWYNQvw8xMQjEgadcYSu2naFOjQgWbYmqOhGa29e3siOZmKPLEvKvTE\nIg3NsM3L0+Pvfwe+/RbQ02RN6HTA3r3AuXPGf1/NmtEvidgftdETiyiVCmi1yQZtzo89loQpU4bg\n+nVg2jTgwgXetj90KDBkCNDAStYAnLdj15JcFy7wWaw7d/LNXjp1Anr1UqCqKhnnz9/9ffEZrUPs\n/SMQQoWeWOZOscvISKkzQ3NI7f3//CdQUsIL3rZtgFIJdO7Mi/6wYUCfPnf3MW2oo7Lu44hgaq7q\nav4OZudOfjt3jr/AxcfzyWd8S4YYqFQN/74IsSfqjCUOUVUFHDhwtxiePw8oFLzwr18/F2q17Tp2\nbfXuoKEO1Li4FKxbt6D2Z/n6a37VPnQov/XtC3jRJRRxABpeSZzKQw/x3a3k8vpX+3v32q5j15bv\nDhrqQN23zxOhofyqffhwYOVKICDA7KiEOAwVeiJEUBDf13TKFECh0MHYrpVHj+oRF8fb9uve/P0N\nP/b1vfs16enZBkUeALTaNGRkpNQr9LdvA5cv8zVlLl+ufzt+3HgHamioHvv20VU7kQ46VYlwM2Yo\ncPasYcdux45JmDVrCLp0MSy+J0/WL8w+PneLfmGh8VO6oMATI0cafl11tfEXjsBAICoKeOIJBdau\nTUZJiWEHakrKECryRFLodCXCPahj934YA27cuHtl/tprOly5Uv+4li31ePllw4LetCng4XG/7x6D\nnj2pA5W4ABtuX2gxJ4lhsT179oiOYDEpZ2esfv6srFwmkyUZ7Kcqk81hWVm5YgI+gKv9/qVG6vlN\nrZ0WT5j64osvEBISAk9PTxw9erTB4zQaDXr06IHw8HBkZGRY+nBOTa1Wi45gMSlnB+rnj4+PwYoV\ncYiLS8GAAamIi0vBihXOexXuar9/qZF6flNZ3HQTFhaGLVu2YNq0aQ0eo9frMWnSJHz11VcIDAxE\n7969MXjwYHTr1s3ShyXkgeLjY5y2sBMigsWFvmvXrg885vDhw+jcuTOCg4MBAGPHjsW2bduo0BNC\niCNZ20Ykl8vZ999/b/RzX3zxBZsyZUrtxx9//DF77bXX6h0HgG50oxvd6GbBzRT3vaKPjY3FpUuX\n6t2/aNEijBgx4n5fCoDP2jIFo1mxhBBiN/ct9DnGZrGYITAwEMXFxbUfFxcXIygoyKrvSQghxDw2\nWaa4oSvyXr16obCwED///DOqqqrw+eefIyEhwRYPSQghxEQWF/otW7agQ4cOOHjwIOLj4zF06FAA\nwIULFxAfHw8A8PLywrp16zB69Gj07NkTkyZNoo5YQghxMOGrV2o0GsycORM6nQ5Tp05FYmKiyDhm\nmTRpElQqFdq2bYvjx4+LjmO24uJivPTSS7h8+TLatGmDiRMnYuLEiaJjmaSiogIDBgxAZWUlfH19\n8ec//xmzZs0SHctser0evXr1QlBQELZv3y46jlmCg4PRrFkzeHp6wtvbG4cPHxYdySw3b97E9OnT\nUVBQgMrKSqxbtw59+/YVHcskp0+fxtixY2s/Pnv2LBYsWAClUmn8C8wfZ2M7Op2OyWQyVlRUxKqq\nqlhERAT78ccfRUYyi0ajYUePHmWhoaGio1jk4sWLLC8vjzHG2K+//sr8/f0l9fu/efMmY4yxiooK\nFhISwgoLCwUnMt+7777Lxo8fz0aMGCE6itmCg4PZ1atXRcew2EsvvcQ++ugjxhhj1dXV7Pr164IT\nWUav17OAgAB27ty5Bo8RupVg3XH23t7etePspaJ///5o2bKl6BgWCwgIQGRkJACgdevW6N27Ny5c\nuCA4len8/thotby8HDqdDj4+PoITmaekpAQ7duzAlClTJDvyTKq5f//9d+zduxeTJk0CwJuZmzdv\nLjiVZb766ivIZDJ06NChwWOEFvrz588bhAsKCsL58+cFJnJfZ86cwYkTJyTz1hUAampqEBERAX9/\nf7z22mv3PdGd0axZs7BkyRI0aiTNrZs9PDwwaNAgREVFYe3ataLjmKWoqKi2uTI0NBRTp07F7du3\nRceyyH/+8x+MHz/+vscIPcNMHWdP7Ku8vBxjx47FsmXL0KRJE9FxTNaoUSPk5+fjzJkzeP/995GX\nlyc6ksmysrLQtm1bREVFSfaqeP/+/cjPz8enn36KRYsWYe/evaIjmUyn0+HIkSMYM2YMjhw5gsrK\nSnzxxReiY5mtqqoK27dvx7PPPnvf44QWehpnL151dTXGjBmDF154ASNHjhQdxyLBwcEYNmwYcnNz\nRUcx2YEDB5CZmYlOnTph3Lhx+Oabb/DSSy+JjmWWdnwzXHTr1g2jR4+WVGdsUFAQWrVqhREjRqBx\n48YYN24cdu7cKTqW2Xbu3ImePXuiTZs29z1OaKGncfZiMcYwefJkhISEYObMmaLjmOXKlSu4fv06\nAODq1avYuXMnwsLCBKcy3aJFi1BcXIyioiL85z//waBBg7Bp0ybRsUx269YtlJWVAQB+/fVX7Nix\nQ1K//4CAAHTu3BmHDh1CTU0NVCoVBg8eLDqW2T777DOMGzfuwQc6qGO4QWq1mkVGRrLQ0FC2YsUK\n0XHMMnbsWNauXTv20EMPsaCgILZu3TrRkcyyd+9e5uHhwSIiIlhkZCSLjIxkO3fuFB3LJAUFBSwq\nKoqFh4czhULB/vWvf4mOZDG1Wi25UTdnz55lERERLCIigg0aNIitXr1adCSznT59mkVHRzOZTMZG\njRrFysvLRUcyS3l5OWvVqhW7cePGA48VPo6eEEKIfUmzu58QQojJqNATQoiLo0JPCCEujgo9IYS4\nOCr0hBDi4qjQE0KIi/v/qZVoR7e5egQAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again we get good agreement if we shift by the difference at the left-most point:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u2 = u + v[0]-u[0]\n",
      "norm(u2 - v, inf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "1.5679084253150677e-08"
       ]
      }
     ],
     "prompt_number": 10
    }
   ],
   "metadata": {}
  }
 ]
}