{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<h1>Fast Fourier Transform</h1>\n",
        "\n",
        "The discrete Fourier transformation interpolates a set of values in the frequency basis given by trigonometric functions. \n",
        "\n",
        "<h2>DFT Matrix</h2>\n",
        "\n",
        "Definte the $n$th root of unity as $\\omega_{(n)} = e^{-2\\pi i/n}$ so that ${\\omega_{(n)}}^n=1$. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "(-1-1.22464679915e-16j) (-0.5-0.866025403784j) (6.12323399574e-17-1j)\n",
            "(1+2.44929359829e-16j) (1+6.10622663544e-16j) (1+2.44929359829e-16j)\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "\n",
        "def omega(n):\n",
        "    return np.exp(-2.*np.pi*1.j/n)\n",
        "\n",
        "print(omega(2),omega(3),omega(4))\n",
        "print(omega(2)**2,omega(3)**3,omega(4)**4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets plot $\\omega_{(2)},\\omega_{(3)},\\ldots,\\omega_{(n)}$ on the complex plane, highlighting in red $\\omega_{(2)},\\omega_{(4)},\\omega_{(8)},\\ldots,\\omega_{(n)}$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ -1.00000000e+00 -1.22464680e-16j  -5.00000000e-01 -8.66025404e-01j\n",
            "   6.12323400e-17 -1.00000000e+00j   3.09016994e-01 -9.51056516e-01j\n",
            "   5.00000000e-01 -8.66025404e-01j   6.23489802e-01 -7.81831482e-01j\n",
            "   7.07106781e-01 -7.07106781e-01j   7.66044443e-01 -6.42787610e-01j\n",
            "   8.09016994e-01 -5.87785252e-01j   8.41253533e-01 -5.40640817e-01j\n",
            "   8.66025404e-01 -5.00000000e-01j   8.85456026e-01 -4.64723172e-01j\n",
            "   9.00968868e-01 -4.33883739e-01j   9.13545458e-01 -4.06736643e-01j\n",
            "   9.23879533e-01 -3.82683432e-01j]\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEACAYAAACUHkKwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXe9vHvL2QhaMAghCQoiTsqKCQYUMdJJCgoEJBx\nHxwQEBDXZxTlUZ8XkEfRmWuUUUSYV8AdnRGERMXd1hccjYQdIsGR1UQIi2wSsp33j+7EiNmru093\n5/e5Li463UXVnYbcnDpVXSXGGJRSqrnCbAdQSgU3LRGllCNaIkopR7RElFKOaIkopRzRElFKOeKV\nEhGRuSKyS0TW1vF6uoj8JCIrPb8e8cZ2lVL2hXtpPfOBZ4GX61nmC2NMlpe2p5QKEF4ZiRhjlgH7\nG1hMvLEtpVRg8eecSB8RWSUi74rIeX7crlLKh7y1O9OQPCDJGPOziFwFLAbO9tO2lVI+5JcSMcYc\nrvF4qYjMEpH2xph9xy8rIvphHqUsMcY0edrBm7szQh3zHiLSqcbjNEBqK5AqxpiA+jV58mTrGTRT\n6GQK1FzN5ZWRiIi8DmQAJ4vIdmAyEOnuA/MP4FoRuR0oA44CN3hju0op+7xSIsaYmxt4/TngOW9s\nSykVWPSM1UbIyMiwHeE3NFPjBGImCNxczSFO9oV8QURMoGVSqiUQEYzliVWlVAukJaKUckRLRCnl\niJaIUsoRLRGllCNaIkopR7RElFKOaIkopRzRElFKOaIlopRyREtEKeWIlohSyhEtEaWUI1oiSilH\ntESUUo5oiSilHNESUUo5oiWilHJES0Qp5YiWiFLKES0RpZQjWiJKKUe0RJRSjmiJKKUc0RJRSjmi\nJaKUcsQrJSIic0Vkl4isrWeZZ0Rks4isFpEe3tiuUso+b41E5gP963pRRK4CzjDGnAWMA2Z7abtK\nKcu8UiLGmGXA/noWGQK87Fn2a6CdiHTyxrZVYNiTX8zGl75hT36x7SjKz/w1J9IZ2FHj6x88z6kQ\nsPyuBbQ5L4nOI6+gzXlJLL9rge1Iyo/CbQdQge/QoUMUFRVRWFj4q993797NweKD8H4OhkoMRwkH\nZOYtPLXlLeKT4klMTCQhIeFXv5988smEhemcfqjwV4n8AJxa4+tTPM/VasqUKdWPMzIyyMjI8FUu\nVUNpaSkbNmwgLy+v+ld+fj6VlZW/KYKEhAS6d+9O8Vc7ieMj2vEzABXAQSLZ3bkH0V1PoqioiC++\n+OJX5XPw4EGSkpJISUkhNTWV1NRUUlJSiI2NrTVXcTFs3QrJydCxo9/ejpDncrlwuVyO1yPGGOdp\nABFJBnKMMd1ree1q4A5jzEAR6QPMMMb0qWM9xluZVP3279/P+++/z+eff05eXh4bNmzgtNNOq/7B\nTk1NpVu3brRr1w4RqXUde/KLaXNeEm04Wv3cz0Tz88ZtdDi39p/4Y8eO8Z///OdXZbV69Wri4uJI\nTU2lT58+DBw4kHPOOYcFC2D0aIiMhNJSmDsXbrrJJ29HiyciGGNq/4uu78954wdWRF4HMoCTgV3A\nZCASMMaYf3iWmQkMAI4AtxpjVtaxLi0RH/r+++/Jzs4mOzubFStWkJGRQWZmJr169aJHjx6ccMIJ\nTV7n8rsW0HPmaMqIIIIyVt05l0ufbdpPekVFBQUFBeTl5bFs2TJycnKIjj6RbduyKC/PAi4GwomO\nhm3bdETiC1ZLxJu0RLxv3bp1LFiwgOzsbIqLixk8eDBZWVn069ePNm3aeGUbe/KL2Z27lbi05DpH\nIE1hjOGVV1Yydmw2x45l456XH0h09DA++WQgF1+s03nepiWifqW0tJSFCxcya9Ysvv/+e/70pz8x\nZMgQ0tLSgmZSs7gYkpLg6FGA7UAOYWGvEx+/nfHjxzJmzBgSEhIspwwdzS0RjDEB9csdSTXXtm3b\nzEMPPWQ6depkMjMzzcKFC01paantWM32+uvGREcb07at+/fXXzdm9erVZty4ceakk04y119/vXG5\nXKaystJ21KDn+dlr8s+sjkRCxIoVK5g2bRrLli1j+PDh3H777XTt2tV2LK+o6+jMgQMHeOWVV5g1\naxYiwv33388tt9xCeLju6jSH7s60IDV/qPbvL+CRRx5h+fLlPPTQQ4wcObJZk6PBzBiDy+ViypQp\nFBcX89hjjzF06FD27BE9NNwEujvTQlQN72NifjCtWo01J554snn88cfNkSNHbEezrrKy0rz33nvm\nwgsvNGee2dtERn5m2rX7ZTdI1Q/dnQl9xcXQpctPlJQ8AfxfYDStW09i+/b2+j9tDbt2VXLqqW9Q\nVvYIcA4wnejoHnpouAHNHYkExzS9AuCVV3I4dux8oBhYA/yFyMj2bN1qN1eg2b49jDZtbga+BQYB\nA6iomMSmTSWWk4UmLZEgsG/fPm655RaeeeZeIiJeB+bi/uQAlJW59/nVL5KT3We3us93vANYS3n5\nd4wenUJubq7VbKFISyTA5eTk0L17d2JjY9mwYS0vvphOdDS0bQvR0e7TwHWI/msdO7rfl1/epzhe\neeVfPProZLKyspg0aRIlJToq8RadEwlQ+/bt45577uHLL79k3rx5pKenV7+mH0hrnNrep927dzNh\nwgQ2btzI/Pnz6d27t76fHnp0JoSsWbPGJCcnmzvvvNMcPnzYdpyQU1lZad544w0TFxdnRo6caVq3\nrtSjOEaPzoSMRYsWMW7cOJ555hlu0o+r+lRu7n/o02cIxlwKPAtEtugP+OnRmSBXWVnJ1KlTuffe\ne1m6dKkWiB+InEFMzL+BH4F+QDEREejRribS84MDwJEjRxgxYgQ//PADubm5xMfH247UIiQnQ1lZ\nDPA28H+ANI4dW0xy8oV2gwUZHYlYtnPnTi699FJiYmJwuVxaIH70y1GcMNq2/V8iIp4gMrIfy5a9\nbTtaUNE5EYu2bNlCZmYm48aN44EHHqjz6mHKt2oendm+PY9Bgwbx17/+lQGp/b16jZRApx/ACzKb\nN2+mX79+PPDAA9xxxx2246gaNm7cSHra73j0yGFupg0RlDbram3BRkskiGzevJnLL7+cqVOnMnr0\naNtx1HH25Bfzw3mnMphjPAKMpeHrxoYCPToTJKp2YbRAAtfu3K0k05pPgUdx33WtjAh25261GyxA\n6dEZP9qxYweZmZk8+OCDWiABLC4tmQhKORP4CMgE4Ch9zkzmm2/0zNbj6UjETw4fPsygQYMYP368\nzoEEuA7ndmTVnXP5mWgSacvbRHFnRDTd+xZwxRXu674u0Jv8VdM5ET+orKzk2muvJTY2lhdeeEGP\nwgSJqivYh5+ZTPe+KygtHQ18BXQJyTNbdU4kgD366KPs2rWr+lqgKjh0OLcj5424iAORHYmOvgr4\nM+570x/RM1tr0BLxsYULFzJv3jwWLlxIVFSU7TiqGX65Psl9QHfgVo4dM+zf7z7HpKXTEvGhNWvW\nMH78eBYvXqxnogaxX85sFWJi/oHINsrLH+f663V+BHROxGf27dtHSkoKTz75JDfccIPtOMoLioth\n1SoYMqSQkpLewCxgcMjMj+icSIC5++67GTJkiBZICOnYEWJjISoqEXgD92loe1r8/IieJ+IDS5Ys\n4euvv2bNmjW2oygv+2V+5FLgZuAuSksXVM+PBPtopDl0JOJle/fu5fbbb2fevHleu1m2Chw1r98a\nE/O/wErKyxe16PkRr8yJiMgAYAbuUpprjHnyuNdHAH8FdnqemmmMmVfHuoJ6TmT48OF06NCBGTNm\n2I6ifKhqfmTw4C8pLf0DsA7oENTzI82dE3G8OyMiYcBM3GcHFwLfiMgSY8y3xy36hjHmbqfbC2S6\nG9NyVM2PREdfQmmpe7cGFhAW5i6XK6+0ndB/vLE7kwZsNsZsM8aU4Z5xGlLLciF9ltX+/ft1N6aF\n+WV+xL1bA4s4cgSGDm1ZuzXeKJHOwI4aX+/0PHe8YSKyWkT+KSKneGG7AWX69OkMGjSIyy67zHYU\n5SdV8yOtW0cDs4H7gVKOHoXRo1vOiWj+OjqTDbxujCkTkbHAS1R9OLIWU6ZMqX6ckZFBRkaGr/M5\nsnPnTubOncu6detsR1F+dtNNcPLJMGzY5Rw5cg7wD+BOWrVyH/YN5LkRl8uFy+VyvB7HE6si0geY\nYowZ4Pl6Eu77VzxZx/JhwD5jzEl1vB50E6u33XYbHTp0YPr06bajKAuKi91HZo4eXQ1cBWwGTmT2\nbBg3znK4JrB2ZTMRaQVswj2yKAJygZuMMfk1lok3xvzoeXwNMNEYc0kd6wuqEvn222+57LLLKCgo\nIDY21nYcZcmcOTB+PLjPHTkX+B9at4bt2wN7NFKTtTNWjTEVwJ3Ah8AG3Edh8kVkqogM8ix2t4is\nF5FVnmVHOt1uoHj44YeZOHGiFkgLl5ICMTEA04C/A8WUlLjLJdTpZ2ccyM3NZdiwYRQUFOgRmRau\nuBi6dAH3fcLvAKKAp4JqNKKfnbFg2rRpPPLII1ogio4d4eGHq776H2A+LWU0oiORZtqyZQsXXXQR\n27dv1xJRwPGjkVtpwymcRxY/RiWzckfHgB+N6EjEz+bMmcOIESO0QFS1mqORKziDjjzG+/Rj07Ek\n/t8doXv2mY5EmqGkpIQuXbqwfPlyzjrrLNtxVAApLoaUU4vZdCyJDI4yGRhIcNy3RkcifvTWW2/R\ns2dPLRD1Gx07wpSRWyklkgm4L1sEoX3fGi2RZpg1axYTJkywHUMFqCH3JBNJKTfgPmnqeyCCMuLS\nku0G8xEtkSZatWoVO3fuZODAgbajqABVdd8aQzQ3EslMwll159yA3pVxQudEmmjixIm0bt2aadOm\n2Y6iAtye/GK+fPtTbn/2z+ws3BnwtwvRORE/yc7OZujQobZjqCDQ4dyODP7v62lzYhtWr15tO47P\naIk0waZNmzh8+DApKSm2o6ggISJkZWWRnZ1tO4rPaIk0QU5ODoMHDw74YakKLFoiqlp2djZZWVm2\nY6ggc+mll7J161Z27tzZ8MJBSEukkfbs2cOaNWvo27ev7SgqyISHh3P11VeTk5NjO4pPaIk00nvv\nvUdmZiatW7e2HUUFoVDepdESaaQvvviCfv362Y6hglRmZibLli2joqLCdhSv0xJppBUrVtCrVy/b\nMVSQat++PXFxcRQUFNiO4nVaIo1QUlJCQUEBF1xwge0oKoilpqaSl5dnO4bXaYk0wtq1azn77LN1\nPkQ5oiXSguXl5emujHJMS6QFy8vLIzU11XYMFeRSUlJYvXo1lZWVtqN4lZZII6xYsUJLRDnWvn17\nOnToEHKTq1oiDTDG8O2339KtWzfbUVQI6N69Oxs3brQdw6u0RBqwb98+oqOj9VqqyisSExMpKiqy\nHcOrtEQaUFRUREJCgu0YKkRoibRAhYWFJCYm2o6hQkRCQgKFhYW2Y3iVlkgDdCSivElHIi1QUVGR\njkSU1+hIpAUqLCzUkYjyGh2J1EFEBojItyJSICIP1vJ6pIi8ISKbReTfItLFG9v1h927d9OpUyfb\nMVSI6NixI3v37g2pT/M6LhERCQNmAv2B84GbRKTrcYuNBvYZY84CZgB/cbpdfzm05yDF/97Bnvxi\n21FUCAgLCyM8PJyysjLbUbzGGyORNGCzMWabMaYMeAMYctwyQ4CXPI/fAjK9sF2fW37XAswnHxD3\n7FTanJfE8rtC936qyn/Cw8MpLy+3HcNrvFEinYEdNb7e6Xmu1mWMMRXATyLS3gvb9pk9+cX0nDma\nSippx8+04Sg9Z47WEYlyLCIiIqRKJNzSduu9XPqUKVOqH2dkZJCRkeHjOL+1O3crnYnEcLQ6bNX9\nVEP1TmbKP0QkID6E53K5cLlcjtfj+A54ItIHmGKMGeD5ehJgjDFP1lhmqWeZr0WkFVBkjImrY30B\ncQe8PfnFtDkvies4ygSC587uKvDFxMRQWFhITEyM7Si/YvMOeN8AZ4pIkohEAjcCx1+RNgcY4Xl8\nHfCpF7brU1X3UxXC+IlofiY6pO+nqvynvLyc8HBbOwHe5/g7McZUiMidwIe4S2muMSZfRKYC3xhj\n3gHmAq+IyGZgL+6iCXiXPnsTUvA6u8+8hJ/vHMOlWiDKC8rLy2nVqpXtGF7jlTo0xrwPnHPcc5Nr\nPD4GXO+Nbflbp1M7ccKFJ+sIRHnFoUOHiIyMJCIiwnYUr9EzVhsQimcYKnuqPkYRSrdi1RJpQCh+\n1kHZE4ofo9ASaUBCQoKORJTXhOIHOrVEGqC7M8qbdCTSAunujPKmULzIlZZIAzp16sTu3btD6lOX\nyp5QvMiVlkgDIiMjiY+PZ9u2bbajqBDw3XffkZycbDuGV2mJNEJKSkpI3rlM+VdZWRkbNmygR48e\ntqN4lZZII4Tq7Q+Vf23YsIGkpCROPPFE21G8SkukEbRElDeE6u1YtUQaoapEAuHTxSp4aYm0YPHx\n8bRp04YtW7bYjqKCmJZIC6e7NMqJsrIy1q9fT8+ePW1H8TotkUbq3bs3y5cvtx1DBamVK1dy2mmn\nhdykKmiJNNrAgQPJzs7WeRHVLDk5OQwaNMh2DJ/QEmmkCy64gIqKCjZu3Gg7igpC2dnZZGVl2Y7h\nE1oijSQiZGVlkZOTYzuKCjJbtmzhxx9/pHfv3raj+ISWSBNkZWWRnX385WOVql/VrkwoXRKxJi2R\nJkhPT2fjxo3s2rXLdhQVREJ5Vwa0RJokMjKSK6+8knfffdd2FBUkDhw4QG5uLldccYXtKD6jJdJE\nw4YNY8ECvZ2mapyFCxdy+eWXc8IJJ9iO4jNaIk10zTXXsHbtWgoKCmxHUQHOGMNzzz3H+PHjbUfx\nKS2RJoqKimLUqFHMnj3bdhQV4L755hv2799P//79bUfxKce30fS2QLmNZn22bt1Kr1692L59O23a\ntLEdRwWokSNHcv755zNx4kTbURrF5m00W5zk5GQuvvhinRtRddq7dy9Llizh1ltvtR3F57REmmnC\nhAk899xzehq8qtX8+fPJysqiQ4cOtqP4nJZIM/Xv358DBw7w9ddf246iAkxFRQXPP/88EyZMsB3F\nL7REmiksLIw///nPTJ061XYUFWBee+014uPjSUtLsx3FLxxNrIpILPAmkARsBa43xhyoZbkKYA0g\nwDZjzNB61hnwE6tVSktL6dq1K/PmzSMjI8N2HBUAjh07xjnnnMOrr77K7373O9txmsTWxOok4GNj\nzDnAp8B/17HcEWNMijGmZ30FEmwiIyOZNm0akyZN0rkRBcDs2bPp1q1b0BWIE05HIt8C6caYXSIS\nD7iMMV1rWe6QMSamkesMmpEIQGVlJSkpKUyePJlrrrnGdhxl0cGDBznrrLP4+OOP6d69u+04TWZr\nJBJnjNkFYIz5EYirY7koEckVkS9FZIjDbQaUsLAwpk+fzkMPPUR5ebntOMqiv/3tb/Tv3z8oC8SJ\n8IYWEJGPgE41nwIM8Egti9c1hEgyxhSJyGnApyKy1hhT51WPp0yZUv04IyMj4OcbBgwYwBNPPMHL\nL7/MqFGjbMdRFuzevZuZM2eyYsUK21EazeVy4XK5HK/H6e5MPpBRY3fmM2PMuQ38mflAjjFmUR2v\nB9XuTJWvvvqKP/zhD6xfv57Y2FjbcZSfjRo1irZt2zJjxgzbUZrN1u5MNjDS83gEsOT4BUTkJBGJ\n9DzuAFwChNw1Bvv06cM111zDvffeazuK8rOlS5fy2WefMW3aNNtRrHA6EmkP/BM4FdiG+xDvTyKS\nCowzxowVkYuBOUAF7tJ62hjzYj3rDMqRCMDhw4e58MILmTFjBoMHD7YdR/nBTz/9RPfu3XnppZfo\n27ev7TiONHckoh/A8zKXy8Xw4cNZt26d7ta0AKNGjSIqKornn3/edhTHtEQCyF133cXBgwd56aWX\nbEdRPrR06VImTJjA2rVriYlp1BkMAU1LJIBU7db8/e9/D9l7jbR0obQbU0VLJMB8/vnn3HTTTeTm\n5nLKKafYjqO8yBjDDTfcQIcOHZg1a5btOF6j1xMJMOnp6dxzzz0MHTqUo0eP2o6jvOixxx5j+/bt\nPPXUU7ajBAQdifiQMYbhw4djjOG1115DpMklrwLM22+/zd13301ubi4JCQm243iVjkQCkIjwwgsv\nUFBQwJNPPmk7jnJo3bp1jB07lkWLFoVcgTjR4Gnvypno6GgWL15M79696datm060Bqk9e/YwZMgQ\nZsyYwUUXXWQ7TkDR3Rk/+eqrrxg8eDAul4vzzz/fdhzVBKWlpfTv35+0tLSQHlHq0Zkg8NprrzFp\n0iRcLhdnnHGG7TiqEcrLy7nxxhspLy9n4cKFIXs/XWh+iejujB/98Y9/5PDhw2RmZvLFF1/QpUsX\n25FUPSoqKhg5ciRHjhxh8eLFIV0gTmiJ+Nm4ceMoKSmhb9++fPbZZ5x66qm2I6laVFRUcNttt1FY\nWMi7775LVFSU7UgBS0vEgnvuuYfKykrS09P55JNPOO2002xHUjWUl5czcuRICgsLyc7OJjo62nak\ngKYlYsl//dd/ERUVRXp6Oh9//DFnn3227UgKKCsr4+abb+bQoUO8++67WiCNoCVi0YQJE2jdujW/\n//3vefPNN0lPT7cdqUXbu3cv1113HTExMSxevJjWrVvbjhQU9GQzy0aNGsXLL7/M9ddfz5w5c2zH\nabHWr19PWloaF110EYsWLdICaQI9xBsgNm/eTFZWFn379mXGjBlERETYjtRiLFmyhDFjxvD0008z\nfPhw23Gs0fNEQsCBAwe4+eab+fnnn/nXv/7VIu7japMxhscff5znn3+eRYsWtZg71tVFPzsTAtq1\na0d2dja9e/cmLS2NlStX2o4Usg4cOMCNN97IkiVLyM3NbfEF4oSWSIBp1aoVTzzxBNOnT2fAgAFM\nnjyZ0tJS27FCygcffMAFF1xAbGwsn3/+OYmJibYjBTXdnQlghYWFjBs3ju3bt/Piiy/Ss2dP25GC\n2oEDB7j//vv56KOPeOGFF+jXr5/tSAFFd2dCUGJiItnZ2dx33330799fRyUOfPDBB3Tv3p1WrVqx\nbt06LRAv0pFIkCgsLGTs2LHs2LGDOXPm0KdPH9uRgkJxcTEPPfSQjj4aQUciIS4xMZGcnBwmTpzI\ntddey7Bhw8jPz7cdK2AdOnSIRx99lHPPPZfo6GjWrl2rBeIjWiJBREQYPnw4mzdv5uKLL+b3v/89\nY8aMYceOHbajBYzS0lKeffZZzjrrLDZt2kRubi7PPPMMbdu2tR0tZGmJBKHo6GgmTpxIQUEBHTt2\npEePHkycOJF9+/bZjmZNZWUlr776Kl27dmXp0qW8//77vPbaa5x++um2o4U8LZEgFhsby/Tp01m3\nbh2HDh3izDPPZMKECaxfv952NL/Zv38/Tz/9NF27duW5555j/vz5vPfee/To0cN2tBZDSyQEJCYm\nMnv2bNatW0dcXBxXXnkl6enpvPnmmyF7NGflypWMHj2a008/nby8PF588UW+/PJL/RCjBXp0JgSV\nlZWxZMkSZs2aRX5+PmPGjGHUqFFBf92SQ4cO8fbbbzNr1iyKiooYP348o0ePJi4uzna0kGDlszMi\nci0wBTgXuMgYU+t52iIyAJiBe+Qz1xhT59VutUS8a+PGjTz//PO8+eabxMfHk5WVRVZWFr169SIs\nLPAHojt37iQnJ4fs7GyWL1/OZZddxvjx47n66qv1coVeZqtEzgEqgTnA/bWViIiEAQVAJlAIfAPc\naIz5to51aon4QEVFBV9//XX1D+S+ffsYNGgQWVlZXH755Zx44om2IwLunKtXr67OuW3bNq6++mqy\nsrLo37+/HmXxIauf4hWRz4D76iiRPsBkY8xVnq8nAaau0YiWiH9899131T+oubm5JCUl0atXL1JT\nU0lNTaVHjx4+L5aKigo2bdrEihUryMvLIy8vjzVr1tC5c+fqgrvkkksID9drZ/lDIJfIH4D+xpix\nnq+HA2nGmLvrWJeWiJ+VlZWxYcMG8vLyqn+gN2zYQFJSEt26daNz584kJCSQkJBAYmJi9e/t2rWr\n99agJSUlFBUVUVRURGFh4a8ef/fdd6xZs4ZOnTr9qrxSUlI46aST/Pjdqyo+KxER+QjoVPMpwAAP\nG2NyPMtoiYSYqmLJz8+vtQSKioooKSkhKiqK8PBwwsPDCQsLo7y8nLKyMsrKyqioqCA+Pr66eGqW\nUHJyshZGgPHZfWeMMVc0L1K1H4CaN1g5xfNcnaZMmVL9OCMjg4yMDIcRVFNFRETQo0ePes+3KCkp\nobS0tLo4KisriYiIqC6VE044QW9iHsBcLhcul8vxery5O3O/MSavltdaAZtwT6wWAbnATcaYWj/4\noSMRpeyw8gE8ERkqIjuAPsA7IrLU83yCiLwDYIypAO4EPgQ2AG/UVSBKqeCjJ5sppQC9FIBSyhIt\nEaWUI1oiSilHtESUUo5oiSilHNESUUo5oiWilHJES0Qp5YiWiFLKES0RpZQjWiJKKUe0RJRSjmiJ\nKKUc0RJRSjmiJaKUckRLRCnliJaIUsoRLRGllCNaIkopR7RElFKOaIkopRzRElFKOaIlopRyREtE\nKeWIlohSyhEtEaWUI1oiSilHtESUUo44KhERuVZE1otIhYik1LPcVhFZIyKrRCTXyTaVUoHF6Uhk\nHXAN8HkDy1UCGcaYnsaYNIfb9DuXy2U7wm9opsYJxEwQuLmaw1GJGGM2GWM2A9LAouJ0WzYF4l+4\nZmqcQMwEgZurOfz1g22AD0TkGxG5zU/bVEr5QXhDC4jIR0Cnmk/hLoWHjTE5jdzOpcaYIhHpCHwk\nIvnGmGVNj6uUCjRijHG+EpHPgPuMMSsbsexk4JAx5qk6XnceSCnVLMaYhqYmfqPBkUgT1LpxEWkD\nhBljDovICcCVwNS6VtKcb0IpZY/TQ7xDRWQH0Ad4R0SWep5PEJF3PIt1ApaJyCrgKyDHGPOhk+0q\npQKHV3ZnlFItl9XDroF6sloTcg0QkW9FpEBEHvRxplgR+VBENonIByLSro7lKkRkpee9WuyjLPV+\n3yISKSJviMhmEfm3iHTxRY4mZhohIrs9781KERnlh0xzRWSXiKytZ5lnPO/TahHpYTuTiKSLyE81\n3qdHGlySgwOIAAADCklEQVSpMcbaL+Ac4CzgUyClnuW+B2IDKRfuAv4OSAIigNVAVx9mehJ4wPP4\nQeCJOpY76OP3psHvG7gdmOV5fAPwRgBkGgE8469/Q55t/g7oAayt4/WrgHc9j3sDXwVApnQguynr\ntDoSMQF6slojc6UBm40x24wxZcAbwBAfxhoCvOR5/BIwtI7lfD0x3Zjvu2bWt4DMAMgEvn9vfsW4\nT2PYX88iQ4CXPct+DbQTkU71LO+PTNDE9ylYziINxJPVOgM7any90/Ocr8QZY3YBGGN+BOLqWC5K\nRHJF5EsR8UWpNeb7rl7GGFMB/CQi7X2QpSmZAIZ5dhv+KSKn+DBPYx2f+wd8+2+osfp4doffFZHz\nGlrYm4d4axWoJ6t5KZdX1ZOptv3SumbEkzzv1WnApyKy1hizxctRmyoQDttnA68bY8pEZCzukZKv\nR0jBKA/3v6GfReQqYDFwdn1/wOclYoy5wgvrKPL8Xiwib+MevjoqES/k+gGoOWF4iue5Zqsvk2cy\nrJMxZpeIxAO761hH1Xu1RURcQE/AmyXSmO97J3AqUCgirYC2xph9XszQ5EzGmJpD+BeAv/gwT2P9\ngPt9quL435BTxpjDNR4vFZFZItK+vr+/QNqdqfNkNRE50fO46mS19bZzAd8AZ4pIkohEAjfi/t/O\nV7KBkZ7HI4Alxy8gIid5siAiHYBLgI1eztGY7zvHkxHgOtwT1L7UYCZP8VYZgvffl7oIdf8bygb+\nBCAifYCfqnZZbWWqOScjImm4TwOp/z8Af85W1zITPBT3PuFRoAhY6nk+AXjH8/g03LPtq3BfemBS\nIOTyfD0A2ARs9nUuoD3wsWd7HwIneZ5PBf7heXwxsNbzXq0BRvooy2++b9xnIQ/yPI4C/ul5/Ssg\n2Q9/Zw1lehz3fz6rgE+As/2Q6XWgEDgGbAduBcYBY2ssMxP3kaU11HOE0l+ZgDtqvE9fAr0bWqee\nbKaUciSQdmeUUkFIS0Qp5YiWiFLKES0RpZQjWiJKKUe0RJRSjmiJKKUc0RJRSjny/wEojL+g38M+\niQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f053f6d3ba8>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from matplotlib import pyplot as pt\n",
        "n=16\n",
        "omegas = np.zeros(n-1,dtype=complex)\n",
        "expomegas = np.zeros(4,dtype=complex)\n",
        "\n",
        "for i in np.arange(0,n-1):\n",
        "    omegas[i] = omega(i+2)\n",
        "    if 2**i <= n/2:\n",
        "        expomegas[i] = omega(2**(i+1))\n",
        "\n",
        "print(omegas)\n",
        "pt.scatter(np.real(omegas),-np.imag(omegas),color='b')\n",
        "pt.scatter(np.real(expomegas),-np.imag(expomegas),color='r')\n",
        "\n",
        "circle= pt.Circle((0,0), radius= 1,fill=False)\n",
        "\n",
        "ax=pt.gca()\n",
        "ax.add_patch(circle)\n",
        "pt.axis('scaled')\n",
        "pt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The **discrete Fourier transform** is defined by a matrix-vector product with a **DFT matrix** $\\boldsymbol D \\in \\mathbb{R}^{n\\times n}$ with $d_{jk}={\\omega_{(n)}}^{jk}$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.+0.j,  1.+0.j,  1.+0.j,  1.+0.j],\n",
              "       [ 1.+0.j,  0.-1.j, -1.+0.j,  0.+1.j],\n",
              "       [ 1.+0.j, -1.+0.j,  1.+0.j, -1.+0.j],\n",
              "       [ 1.+0.j,  0.+1.j, -1.+0.j,  0.-1.j]])"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def DFT(n):\n",
        "    v = np.linspace(0,n-1,n) # v=[0,1,...,n-1]\n",
        "    return omega(n)**np.outer(v,v)\n",
        "\n",
        "np.asarray(np.real(DFT(4)),dtype=int)+1j*np.asarray(np.imag(DFT(4)),dtype=int)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The DFT matrix is symmetric (but not Hermitian) and with appropriate scaling, unitary i.e.  $n\\boldsymbol D^{-1} = \\boldsymbol D^{H}$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "2.2721682842512306e-12"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import numpy.linalg as la\n",
        "n=27\n",
        "la.norm(n*np.eye(n) - DFT(n) @ DFT(n).conj())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<h2>FFT Algorithm</h2>\n",
        "\n",
        "The **fast Fourier transform** exploits the special structure of the DFT matrix, computing the disrete Fourier transform of $\\boldsymbol v$, given by $\\boldsymbol u = \\boldsymbol D \\boldsymbol v$, in $O(n\\log n)$ operations."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "3.8305731252303524e-14"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def fft(v):\n",
        "    n = v.size\n",
        "    if n is 1:\n",
        "        return v\n",
        "    u = fft(v[::2]) # compute FFT of [v_0,v_2,...v_{n-1}]\n",
        "    w = fft(v[1::2]) # compute FFT of [v_1,v_3,...v_n]\n",
        "\n",
        "    # scale w by twiddle factors [omega_(n)^0,...omega_(n)^(n/2-1)]\n",
        "    t = np.asarray([omega(n)**i for i in range(n//2)]) \n",
        "    z = w*t\n",
        "\n",
        "    return np.concatenate([u+z,u-z])\n",
        "\n",
        "n = 16\n",
        "v = np.random.random(n)+1.j*np.random.random(n)\n",
        "la.norm(fft(v)-DFT(n)@v)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Two recursive calls on vectors of half the size are made, and $O(n)$ additions/products are done to apply the twiddle factors and add $z$ to $u$, for a cost of\n",
        "$$T(n)=2T(n/2)+O(n)=O(n \\log n).$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true
      },
      "source": [
        "<h2>Application to the 1D Poisson Equation</h2>\n",
        "\n",
        "The solutions to the 1D Poisson equation are waves that correspond to Fourier modes, and the Fourier transform provides solutions to it directly when applied to the continuous equation or in the discretized form. From the linear algebra perspective, the DFT matrix is related to the eigenvectors of the tridiagonal differencing matrix (also linear stiffness matrix) $\\boldsymbol T$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 2., -1.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [-1.,  2., -1.,  0.,  0.,  0.,  0.],\n",
              "       [ 0., -1.,  2., -1.,  0.,  0.,  0.],\n",
              "       [ 0.,  0., -1.,  2., -1.,  0.,  0.],\n",
              "       [ 0.,  0.,  0., -1.,  2., -1.,  0.],\n",
              "       [ 0.,  0.,  0.,  0., -1.,  2., -1.],\n",
              "       [ 0.,  0.,  0.,  0.,  0., -1.,  2.]])"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n=7\n",
        "T = 2*np.eye(n)-np.diag(np.ones(n-1),1)-np.diag(np.ones(n-1),-1)\n",
        "T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The eigenvectors of $\\boldsymbol T$ are $\\boldsymbol X$ with eigenvalues $\\boldsymbol d$ where\n",
        "$$x_{jk} = \\sqrt{\\frac{2}{n+1}}\\sin(jk\\pi/(n+1)), \\quad d_{k} = 2(1-\\cos(\\pi k/(n+1)).$$\n",
        "The eigenvector matrix $\\boldsymbol X$ is the imaginary part of a submatrix of the $2(n+1)$-dimensional DFT matrix. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "2.4218421576539674e-15"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "X = - np.sqrt(2/(n+1.)) * np.imag(DFT(2*(n+1))[1:n+1,1:n+1])\n",
        "d = 2.*(1.-np.cos(np.pi*np.arange(1,n+1)/(n+1)))\n",
        "\n",
        "#check that we indeed have the eigenvalue decomposition of T\n",
        "la.norm(T@X-X@np.diag(d))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Thus we can apply $\\boldsymbol X$, which is called the **discrete sine transform** via FFT."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.0024720669331688e-15"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def fast_DST(v):\n",
        "    # setup a padded vector [0,v,0,...,0] of dimensions 2(n+1)\n",
        "    w = np.concatenate([np.asarray([0.j]),v,0.j*np.zeros(v.size+1)])\n",
        "\n",
        "    u = fft(w) # compute FFT of padded vector\n",
        "    z = u-u.conj() # extract only imaginary part\n",
        "\n",
        "    # return rescaled subvector\n",
        "    return (-1.j*np.sqrt(1./(2.*(v.size+1.))))*z[1:1+v.size]\n",
        "\n",
        "n = 7 # so that DFT dimension is 2(n+1)=16\n",
        "X = - np.sqrt(2./(n+1.))*np.sin(np.pi*np.outer(np.arange(1,n+1),np.arange(1,n+1))/(n+1.))\n",
        "\n",
        "v = np.random.random(n)\n",
        "\n",
        "la.norm( X @ v - fast_DST(v))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true
      },
      "source": [
        "Now, to solve the linear system $$\\boldsymbol T \\boldsymbol x = \\boldsymbol b$$ it suffices to compute $$\\boldsymbol x = \\boldsymbol X \\text{diag}(\\boldsymbol d)^{-1} \\boldsymbol X^{-1} \\boldsymbol b,$$\n",
        "and since the sine transform is orthogonal, we simply have $$\\boldsymbol x = \\boldsymbol X \\text{diag}(\\boldsymbol d)^{-1} \\boldsymbol X^{T} \\boldsymbol b.$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "3.820199830714189e-15"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n = 7\n",
        "b = np.random.random(n)\n",
        "T = 2*np.eye(n)-np.diag(np.ones(n-1),1)-np.diag(np.ones(n-1),-1)\n",
        "x = la.solve(T,b)\n",
        "\n",
        "d = 2.*(1.-np.cos(np.pi*np.arange(1,n+1)/(n+1)))\n",
        "X = - np.sqrt(2./(n+1.))*np.sin(np.pi*np.outer(np.arange(1,n+1),np.arange(1,n+1))/(n+1.))\n",
        "\n",
        "la.norm(x - fast_DST(np.diag(1./d) @ X.T @ b))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<h2>Application to the 2D Poisson Equation</h2>\n",
        "\n",
        "Lets now consider the discretized 2D Poisson PDE, where the linear system takes on a Kronecker-product form."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<matplotlib.image.AxesImage at 0x7f053eda1518>"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD7CAYAAABKWyniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADW1JREFUeJzt3V+opPV9x/H3Z7EVIUE3hlXsJjGlGENIMTemwQvHlDbR\nQIwXDQkIGrxsiZASskkv9njRUgsShEJyEROWQoythrgtFI3oUHqhDXE3SvyzIRe2pt1TKG6DBMRk\nv704j9vj2Tl75syZeWae83u/YGDm2Zl5vnOYz36f329+80yqCkltObDsAiT1z+BLDTL4UoMMvtQg\ngy81yOBLDeol+Ek+meSlJKeSfKWPfe5WkgeSrCd5btO2g0keT/JykseSXLrMGidJcjjJk0l+muT5\nJF/stq9s7UkuTvJMkhNdzUe77Vcnebp7nzyY5KJl17pVkgNJnk1yvLu98jVPsvDgJzkA/C3wCeBD\nwOeTXLvo/c7gO2zUuNkR4Imq+gDwJPDV3qva2a+BL1XVh4CPAX/a/X1XtvaqegO4qao+AlwH3Jzk\no8C9wH1VdQ1wBrhriWVu527ghU23h1Dzefro+NcDP6uqV6rqTeB7wK097HdXqupfgde2bL4VONZd\nPwZ8pteiplBVp6vqZHf9deBF4DArXntV/aq7ejFwEVDATcAj3fZjwG1LKG1bSQ4DtwDf2rT546xw\nzdvpI/i/A/zHptuvdtuG4FBVrcNGwIBDS67ngpJczUYHfRq4YpVr7w6ZTwCngR8CPwfOVNXZ7i6v\nAlctq75tfB34Mhv/SZHkcuC1Fa95Iif3dmdl1zcneQfwMHB31/m31rpStVfV2e5Q/zAbR4WrOPw7\nJ8mngPXu6Cqb/2lJJe1JHxMRvwDeu+n24W7bEKwnuaKq1pNcCfz3sguapJtQehj4u6p6tNs8iNqr\n6pdJxmzMT1yW5EDXQVftfXID8OkktwCXAO8E7gcuXeGat9VHx/8R8HtJ3pfkt4HPAcd72O8swtv/\nBz8O3NldvwN4dOsDVsS3gReq6v5N21a29iTvfutThiSXAH/ExoTZU8CfdHdbqZqr6mtV9d6q+l02\n3sNPVtXtrHDNF1RVC78AnwReBn4GHOljnzPU+F3gP4E3gH8HvgAcBJ7oan8cuGzZdU6o+wbgN8BJ\n4ATwbPf3fteq1g58uKvzJPAc8Bfd9vcDzwCngIeA31p2rdvUfyNwfEg1b72kK15SQ5zckxpk8KUG\n7Sn4Q1iKK+l8M4/xu6W4p4A/ZGNS7EfA56rqpfmVJ2kR9tLxB7EUV9L59rKAZ9JS3Ou33imJHxtI\nS1JVE1cWrsTk3oI/c93T/o4ePbr0z1x3exlizUOte5VrvpC9dPwhL8WV9p3xeMx4PJ7qvnsJ/rml\nuMB/sbGM8fN7eD5JezAajRiNRudu33PPPdved+bgV9VvkvwZG8tBDwAPVNWLszxXcv4wZKdDlWlN\nep6t+7vQvjb/IYdiiDXDMOseYs2wh4/zpt7BjJN7i6xrN8GXhioJtcqTe5L6ZfClBvVyRtCth9KT\nxvRbLfJwfJp6PPzXfmbHlxpk8KUGGXypQQZfatBSfu5nHpN9k55nXvVM2p+TfdpP7PhSgwy+1CCD\nLzVoJX7Sd5ox9iQu8pFmY8eXGmTwpQYZfKlBBl9q0EpM7k3iIh9pcez4UoMMvtQggy81aGXH+FvN\nushnkVzko6Gy40sNMvhSgwy+1CCDLzWol8m9RS10meZ5ppkA7HORzzz3J83Kji81yOBLDTL4UoMM\nvtSglfjtvGVPdvmtPrXGji81yOBLDTL4UoNW9ie0+l5Us9N9PHW39hM7vtQggy81aMfgJ3kgyXqS\n5zZtO5jk8SQvJ3ksyaWLLVPSPE3T8b8DfGLLtiPAE1X1AeBJ4KvzLkzS4mTKb7i9D/jHqvr97vZL\nwI1VtZ7kSmBcVddu89iaZaKqz0Uus57Cq8+anOzTbiWhqia+uWcd4x+qqnWAqjoNHJq1OEn9m9fH\neRdsR2tra+euj0YjRqPRnHYr6S3j8ZjxeDzVfWc91H8RGG061H+qqj64zWM91J+Bh/raq3kc6qe7\nvOU4cGd3/Q7g0Zmr20ZVve2S5LzLovY1bcgWVc+kmhb5+tWeHTt+ku8CI+ByYB04CvwA+AfgPcAr\nwGer6sw2j5+p4094nvO2rdpRQN/1eBSgC7lQx5/qUH+POzf4c2DwtVuLmNWXNGAGX2rQvvrtvL6/\nQbdTPfOsyTP5aJ7s+FKDDL7UIIMvNWgwY/xJhnAmn0XyTD6alR1fapDBlxpk8KUGGXypQYOe3Nuq\n79+nn/IrzXN5nmn0/fo1XHZ8qUEGX2qQwZcaZPClBu2ryb1JhvCttmWvOFz261f/7PhSgwy+1CCD\nLzVo34/xJxnCt/qWfXYhx/37mx1fapDBlxpk8KUGGXypQU1O7m3lqbtd5NMaO77UIIMvNcjgSw1y\njL8NF/m4yGc/s+NLDTL4UoMMvtQggy81yMm9KbnIx0U++4kdX2qQwZcatGPwkxxO8mSSnyZ5PskX\nu+0Hkzye5OUkjyW5dPHlSpqH7DQuS3IlcGVVnUzyDuDHwK3AF4D/qaq/SfIV4GBVHZnw+Gpl7Nf3\nIpc+f55rGi7yWS1JqKqJb5IdO35Vna6qk93114EXgcNshP9Yd7djwGfmU66kRdvVGD/J1cB1wNPA\nFVW1Dhv/OQCH5l2cpMWY+uO87jD/YeDuqno9ydZjuG2P6dbW1s5dH41GjEaj3VUpaUfj8ZjxeDzV\nfXcc4wMkuQj4J+Cfq+r+btuLwKiq1rt5gKeq6oMTHusYv8f99bn/rRzjr5Y9jfE73wZeeCv0nePA\nnd31O4BHZ65wn6iq8y5Jzrsscn/T7L/Peha5P81umln9G4B/AZ5n43C+gK8B/wb8PfAe4BXgs1V1\nZsLjm+n4kyy7Cy77qGDZr79lF+r4Ux3q73HnBn8Lg9/u+6FP8zjUl7SPGHypQX47b8GG8K02v9XX\nHju+1CCDLzXI4EsNcoy/BJ6621N3L5sdX2qQwZcaZPClBhl8qUFO7q0AT93tIp++2fGlBhl8qUEG\nX2qQY/wV5SIfF/kskh1fapDBlxpk8KUGGXypQU7uDYSLfFzkM092fKlBBl9qkMGXGuQYf8CGsMhn\nkVzkMzs7vtQggy81yOBLDTL4UoOc3NtHpp2AW+SE3zT7n+V5ptH36x8yO77UIIMvNcjgSw0y+FKD\nnNzb54bwrbZlrzhc9utfBju+1KAdg5/k4iTPJDmR5PkkR7vtVyd5OsmpJA8m8ehBGogdg19VbwA3\nVdVHgOuAm5N8FLgXuK+qrgHOAHcttFJJczPVoX5V/aq7ejEb8wIF3AQ80m0/Btw29+q0EFX1tkuS\n8y6L2te04+lF1TOppkW+/lU1VfCTHEhyAjgN/BD4OXCmqs52d3kVuGoxJUqat2k7/tnuUP8wcD1w\n7UKrkrRQu5qQq6pfJhkDHwMuS3Kg6/qHgV9s97i1tbVz10ejEaPRaJZaJV3AeDxmPB5Pdd/sNOZK\n8m7gzar63ySXAI8Bfw3cAXy/qh5K8g3gJ1X1zQmPrxY/Jx2Svr/IMssYuu969sN7NglVNfGPPU3w\nP8zG5N2B7vJQVf1lkvcD3wMOAieA26vqzQmPN/gD1Ocil1kn0/qsaYjv4T0Ffw47N/gDZPD3d/Bd\nuSc1yOBLDXKZrSYawqm7l/0TYkM8/H+LHV9qkMGXGmTwpQYZfKlBTu5pKn2fyWaaybWd6plnTfvt\nTD52fKlBBl9qkMGXGuQYXzMbwiKfRRryIh87vtQggy81yOBLDTL4UoOc3NPc9P379NM8zzQTgH1P\nQK7ChJ8dX2qQwZcaZPClBhl8qUFO7mmhhvCttmWvOFzG67fjSw0y+FKDDL7UIMf46t0QvtW37LML\nLXrcb8eXGmTwpQYZfKlBBl9qkJN7WjpP3d3/Ih87vtQggy81yOBLDXKMr5XkIp/Fvn47vtQggy81\naOrgJzmQ5Nkkx7vbVyd5OsmpJA8mcdggDcRuOv7dwAubbt8L3FdV1wBngLvmWZikxZkq+EkOA7cA\n39q0+ePAI931Y8Bt8y1N+n9Vdd4lydsui9zfNLbWM8+a5v36p+34Xwe+DBRAksuB16rqbPfvrwJX\n7WrPkpZmx3F5kk8B61V1Mslo8z9Nu5O1tbVz10ejEaPRaNv7SprNeDxmPB5Pdd/sdBiT5K+A24Ff\nA5cA7wR+APwxcGVVnU3yB8DRqrp5wuNr2SdT1P7U50krZz1s77OmSZ/9V9XEwncM/pYnuhH486r6\ndJKHgO9X1UNJvgH8pKq+OeExBl+96PtMNn3+PNc0tnn9E4vcy+f4R4AvJTkFvAt4YA/PJalHu+r4\nM+3Ajq+e2PH76fiSBsrgSw1yma32jb5/n36a5+lzOLCbMwvZ8aUGGXypQQZfapDBlxrk5J72tVX9\nffrN/O08Sb0w+FKDDL7UIMf4ak7rp+4GO77UJIMvNcjgSw0y+FKDnNxT8/pe5LObb9Fd6D7+dp6k\nXTH4UoMMvtQgx/jSBPt9kY8dX2qQwZcaZPClBhl8qUFO7klTGOoin+3Y8aUGGXypQQZfapBjfGlG\nQ1jks53eOv54PO5rV3M1xLqHWDMMs+4h1gwGf0dDrHuINcMw6x5izeAYX2qSwZcalEX/VE+S5f4+\nkdSwqpo4A7jw4EtaPR7qSw0y+FKDDL7UIIMvNcjgSw36P4UiiKq4dMuUAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f053eda17f0>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "n=7\n",
        "T = 2*np.eye(n)-np.diag(np.ones(n-1),1)-np.diag(np.ones(n-1),-1)\n",
        "A = np.kron(np.eye(n),T)+np.kron(T,np.eye(n))\n",
        "pt.spy(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The 1D finite differencing matrix $\\boldsymbol T$, gives us the following equations for the unknowns $\\boldsymbol U\\in\\mathbb{R}^{m\\times m}$ for an $m\\times m$ spatial mesh,\n",
        "$$\\boldsymbol T \\boldsymbol U + \\boldsymbol U \\boldsymbol T= \\boldsymbol F,$$\n",
        "where we've absorbed the discretization width $h$ into $\\boldsymbol F$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "5.9061360594505111e-15"
            ]
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "F = np.random.random((n,n))\n",
        "U = la.solve(A,F.reshape(n*n)).reshape((n,n))\n",
        "la.norm(T@U+U@T-F)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Using the eigenvalue decomposition of $\\boldsymbol T=\\boldsymbol X \\text{diag}(\\boldsymbol d) \\boldsymbol X^{-1}$, which we know analytically, we can write the frequency space form of the linear system as\n",
        "$$\\text{diag}(\\boldsymbol d) \\boldsymbol X^{-1} \\boldsymbol U \\boldsymbol X + \\boldsymbol X^{-1} \\boldsymbol U \\boldsymbol X\\text{diag}(\\boldsymbol d)  =  \\boldsymbol X^{-1} \\boldsymbol F \\boldsymbol X,$$\n",
        "so we compute the 2D sine transform of $\\boldsymbol F$ and solve a diagonal linear system for $\\boldsymbol V=\\boldsymbol X^{-1} \\boldsymbol U \\boldsymbol X$, which is the 2D sine transform of the solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "7.115828086306871e-15"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "X = - np.sqrt(2/(n+1.)) * np.imag(DFT(2*(n+1))[1:n+1,1:n+1])\n",
        "d = 2.*(1.-np.cos(np.pi*np.arange(1,n+1)/(n+1)))\n",
        "\n",
        "# transform F\n",
        "sF = X.T @ F @ X\n",
        "\n",
        "# perform pointwise product to solve for U\n",
        "sU = sF / np.add.outer(d,d)\n",
        "\n",
        "# backtransform U and compare to naive method\n",
        "la.norm(U - X @ sU @ X.T)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}