{"nbformat_minor": 0, "metadata": {"signature": "sha256:5cebfb4bf1ec0c5df4d88147b485c984f0ff275e4e8a8181703adced7f0cb936", "name": ""}, "worksheets": [{"metadata": {}, "cells": [{"source": ["# Fourier Interpolation"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["import numpy as np\n", "import numpy.linalg as la\n", "import matplotlib.pyplot as pt"], "prompt_number": 2, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["Let's fix the number of points and the number of functions as `n`.\n", "\n", "Make sure `n` is odd."], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["n = 5\n", "\n", "assert n % 2 == 1"], "prompt_number": 49, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [], "metadata": {}, "input": ["x = np.linspace(0, 2*np.pi, n, endpoint=False)"], "prompt_number": 50, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["Next, fix the values of $k$ in $\\cos(kx)$ as `cos_k` and in $sin(kx)$ as `sin_k` so that there are exactly $n$ altogether:"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["cos_k = np.arange(0, n//2 + 1, dtype=np.float64)\n", "sin_k = np.arange(1, n//2 + 1, dtype=np.float64)"], "prompt_number": 51, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [{"text": ["[ 0.  1.  2.]\n", "[ 1.  2.]\n"], "output_type": "stream", "stream": "stdout"}], "metadata": {}, "input": ["print(cos_k)\n", "print(sin_k)"], "prompt_number": 52, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["Next, build the generalized Vandermonde matrix.\n", "\n", "Make sure to order the matrix by increasing $k$:"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [{"prompt_number": 53, "text": ["array([[ 1.        ,  0.        ,  1.        ,  0.        ,  1.        ],\n", "       [ 1.        ,  0.95105652,  0.30901699,  0.58778525, -0.80901699],\n", "       [ 1.        ,  0.58778525, -0.80901699, -0.95105652,  0.30901699],\n", "       [ 1.        , -0.58778525, -0.80901699,  0.95105652,  0.30901699],\n", "       [ 1.        , -0.95105652,  0.30901699, -0.58778525, -0.80901699]])"], "metadata": {}, "output_type": "pyout"}], "metadata": {}, "input": ["V = np.zeros((n,n))\n", "V[:, ::2] = np.cos(cos_k*x[:, np.newaxis])\n", "V[:, 1::2] = np.sin(sin_k*x[:, np.newaxis])\n", "V"], "prompt_number": 53, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["Now let's try and do interpolation with this. Here are a few functions:"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["if 1:\n", "    def f(x):\n", "        return x\n", "elif 0:\n", "    def f(x):\n", "        return np.abs(x-np.pi)\n", "elif 1:\n", "    def f(x):\n", "        return (x<=np.pi).astype(np.int32).astype(np.float64)"], "prompt_number": 54, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["Find the coefficients as `coeffs`:"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["coeffs = la.solve(V, f(x))"], "prompt_number": 55, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [], "metadata": {}, "input": ["plot_x = np.linspace(0, 2*np.pi, 1000)"], "prompt_number": 56, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [], "metadata": {}, "input": ["interpolant = 0 * plot_x\n", "for i, k in enumerate(cos_k):\n", "    interpolant += coeffs[2*i] * np.cos(k * plot_x)\n", "for i, n in enumerate(sin_k):\n", "    interpolant += coeffs[2*i+1] * np.sin(n * plot_x)"], "prompt_number": 57, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [{"prompt_number": 58, "text": ["[<matplotlib.lines.Line2D at 0x7f9157556470>]"], "metadata": {}, "output_type": "pyout"}, {"text": ["<matplotlib.figure.Figure at 0x7f915c0daeb8>"], "metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEACAYAAACqOy3+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0lNX9x/H3TcKigEhFRQUVEEW0GEUBpUIEMVCoCmpV\nCioKWhVwr2sr1lq1am2RRf0huFF6BMSNmgCBYRPZBWUrAawQFdkhQAKZ3N8fN9BAk5DJLM88M5/X\nOTkMM88895sc+ObOvd97r7HWIiIi/pPidQAiIlI1SuAiIj6lBC4i4lNK4CIiPqUELiLiU0rgIiI+\nFVYCN8acY4xZUuprpzFmUKSCExGR8plI1YEbY1KAPKC1tXZDRG4qIiLliuQQypXAWiVvEZHYiGQC\nvwn4RwTvJyIiFYjIEIoxpjpu+KSFtXZz2DcUEZGjSovQfboCi8pK3sYYbbYiIlIF1lpT0euRGkK5\nGRhbQRC+/Xr66ac9jyFZ4/dz7Irf+y+/x18ZYSdwY0wt3ATmh+HeS0REKi/sIRRr7R6gfgRiERGR\nEGgl5lFkZGR4HUJY/By/n2MHxe81v8dfGRFbyFNuA8bYaLchIpJojDHYGE1iiohIjCmBi4j4lBK4\niIhPKYGLiPiUEriIiE8pgYuI+JQSuIiITymBi4j4lBK4iIhPKYGLiPiUEriIiE8pgYuI+JQSuIiI\nTymBi4j4lBK4iIhPKYGLiPiUEriIJJXCwkKmTp3KunXrvA4lbJE41Ph4Y8x4Y8xKY8wKY0zbSAQm\nIhJJ1lqWLl3KsGHDyM/P58QTT/Q6pLCFfaSaMeYdYIa1dpQxJg2oZa3dWep1HakmIp7Ky8sjKyuL\n4uJiunbtSsOGDb0O6agqc6RaWAncGFMXWGKtbVLBNUrgIuKZYDDI6NGjadWqFenp6RhTYU6MG7FI\n4OnAG8AK4AJgEXCftXZvqWuUwEXEU9Za3yTug2JxqHEacBEw3Fp7EbAHeCzMe4qIRJTfkndlpYX5\n/o3ARmvtgpK/j6eMBD548OBDjzMyMsjIyAizWRGRw23bto0vvviCLl26kJYWbmqLvUAgQCAQCOk9\nkZjEnAn0s9b+2xgzGDjGWvtoqdc1hCIiUVNYWMisWbNYvHgxl112GZdeeimpqalehxW2qI+BlzRy\nATASqA6sBfqqCkVEos1ay7Jly8jJyaFJkyZ06tSJOnXqeB1WxMQkgVciCCVwEYm4devWkZOT45uy\nwFApgYtIwjqYVxJ1grIyCdx/I/0iIiRu4g6F9kIRkbiWm5vL4sWLvQ4jLqkHLiJxadu2bWRnZ7Nl\nyxa6dOnidThxSWPgIhJXSpcFtmvXjjZt2viyrjtcmsQUEd+ZOHEixpiEKwsMlRK4iPhOMBhMiIU4\n4YrFXigiIhGl5F15yTewJJKAioshNxe++Qa2boX9++HEE6FxY0hPh2rVvI7wcMFgkPnz59O8eXPq\n1avndTi+pQQu4mNz58Lo0fDxx3DssdCypUvc1avD5s2wejV8+y1ceSX07w9XXQVed3Bzc3PJysqi\nXr16NG/e3NtgfE5j4CI+NHs2PPYYfP89/Pa3cP310KScY1W2boUPP4TXX4cDB+AvfwEvqvIOlgVu\n3ryZLl260KxZMy3GqYAmMUUSzObN8NBDMH06vPAC3HgjVLbCzlr46CP43e+gVSsYMQJiNXpRWFjI\nsGHDaN26NW3btk3KssBQKYGLJJBAAHr3dkn7mWegdu2q3WffPnj0UZfMP/oILrooomGW68CBA1SL\nt8H4OKYELpIArIWXX4a//hXefhsyMyNz3w8/hLvugpEj4ZprInNPiRxtZiXic0VFMGiQG/NesAAi\nuWtqz57QqBFcfTUUFLiefbjy8/NZvnw5bdq0Cf9mclRK4CJx6sAB6NULduyAWbOgbt3It3HJJZCd\n7Xr1KSlwww1Vu08wGGTevHnMnj2b9PR0iouLSUnRMpNo0xCKSBwqKnLJe88eN9RRo0Z021u6FDp3\nhgkT4PLLQ3vvmjVryM7Opl69emRmZlK/fv3oBJlkNAYu4kPFxW6ycvt2mDgRataMTbtTpkCfPjBz\nJpx9duXes2zZMmbMmEFmZqbKAiNMCVzEhx56CBYudEMbsUreB40cCS+95NqvzD5SRUVFACoLjIJY\nHWr8LbALCAIHrLWtj3hdCVykkoYOhWHD4IsvYlejfaQ774T8fBgzBtSh9k6sEvh6oJW1dls5ryuB\ni1TCp5+6sr45c9weJl7Ztw9anzuJi2oNofGJhRTVqMFFvXpxUYcOnHnmmd4FlmRiWUao39MiYVi9\nGu64wyVxL5M3wIJpk+hi7+OlFWsPPXfPokXs+eMfOfOeezyMTI4UiTofC0w1xiw0xvSPwP1Eksqe\nPXDddfCnP0E8lE9PHjKEl75be9hzw7du5d+ffOJRRFKeSPTA21lrfzDGnAhMMcasstbOKn3B4MGD\nDz3OyMggIyMjAs2K+J+1bsz5kkvcboHxIK2wsMznUwsKYhxJcgkEAgQCgZDeE3YCt9b+UPLnZmPM\nRKA1UG4CF5H/Gj4cli93k5bxMmFYVE7ReTDWJTFJ5sjO7TPPPHPU94Q1hGKMOdYYU6fkcS3gKuDr\ncO4pkiy+/hoGD4bx491e3vHiqkGDeLJp08Oe63tsUzoPHOhRRFKecHvgJwMTS4r304Ax1trJYUcl\nkuAKCuA3v4EXX4SzzvImBmst33zzDc2bNz9sl8D23boB8PvXXiO1oIAD1Wvy5YqB9Czu5k2gUi4t\n5BHxwEMPuZNyxo/3ZugkLy+PrKwsiouLueGGGzj++OMrvH7KFDdWv2IFHHNMjIJMclqJKRKHpk6F\n225z+4+ccEJs287PzycnJ4fc3Fw6duxIenp6pZe/33CDO7Lt97+PcpACKIGLxJ0dO+DnP4dRo9zm\nUbG0c+dO3njjDdLT0+nQoQM1Qtwha+1aV+a4YgWcdFKUgpRDlMBF4ky/fu6E+BEjvGl/586d1A1j\nX9r773c7JQ4dGsGgpExK4CJxZOpUuP12+OYbOO44r6Opmi1boHlzmDsXmjXzOprEVpkErh3XRWJg\nzx43CfjGG9FP3oWFheTm5kbl3vXruwnYJ56Iyu0lROqBi8TA/ffDtm3w7rvRa8Nay7Jly8jJyaFZ\ns2Z07949Kvtz790LTZu67W5btoz47aWEhlBE4sAXX8D117uFO9GqOildFti1a1caRvLwzDK88oob\nRhk/PqrNJDUlcBGPFRZCejo8+6xL4tGwcOFCZsyYEXJZYDj27HG98ClTXFWNRJ4SuIjHnn0WFi92\nR6NFS35+PmlpadSM8V4lr7wCX34J48bFtNmkoQQu4qGDddOLF8Ppp3sdTeSpFx5dqkIR8Yi1MGgQ\nPPJI5JL3tm3b2L59e2RuFgG1armKlOef9zqS5KWTSEWi4OOPYf36yAydFBYWMmvWLBYvXkz37t2p\n59VhmWW46y53gtC334JOW4s9DaGIRNiePdCiBbzzDoRzdknpssAmTZrQqVMn6lTmqPgYe/RRt7vi\n3//udSSJRWPgIh547DHIy4P33qv6Pay1vPfeexQWFsakLDAc338P558Pa9bEfnOuRKYELhJjK1ZA\nhw6u5rtBg/DutWnTJk466aSYlAWG6/bboUkTeOopryNJHErgIjFkLVxxhav3HjDA62hia8UK6NjR\njftrv/DIUBWKSAyNGQO7dsHdd4f2vg0bNuD3Tk6LFu5g5mhuFSD/Sz1wkQjYscMlsYkTXe13ZWzb\nto3s7Gy2bNlC3759qV27dnSDjLJp01zp5Ndfx88BzX6mIRSRGBk4EPbvd7sNHk3pssB27drRpk0b\n0tL8X9FrrVvQM2SIG06R8MQsgRtjUoGFwEZr7a+OeE0JXBLaokXQrRssX370KoyffvqJ999/P67L\nAsPx+utul8Jobh2QLGKZwB8EWgF1rLVXH/GaErgkrGAQLr3UjXv37Xv064uKiti0aROnnXZa9IPz\nQH4+nHGG2z7gjDO8jsbfYjKJaYxpCPwSGAlo5EuSysiRUL063Hpr5a5PS0tL2OQNULu2+1kMH+51\nJMkh7B64MWYc8GfgOOBhDaFIsti8Gc47zx2VduTBBsFgkO3bt1O/fn1vgvPQ2rXQti18951KCsNR\nmR54WDMnxpjuwE/W2iXGmIzyrhs8ePChxxkZGWSEs75YJE48+ij07v2/yTs3N5esrCwaNWrENddc\n401wHmra1FXijB3rFvhI5QQCAQKBQEjvCasHboz5M9AHKAJq4nrhE6y1t5S6Rj1wSThz5sCNN7oF\nLAfPuCxdFpiZmUmzZs18sYoyGrKz4fHH3QRvkv4IwhbTMkJjTAc0hCJJoKgILroInnzSJXGA+fPn\nEwgEEqosMBzFxXDWWfDBB3DxxV5H409RH0IpgzK1JLzXXoOTT4Zf//q/z51++uncfffdCVcWWFUp\nKdC/P7z5phJ4NGkhj0glzJw0iclDhhDcVciMhTUYMHQQve7q5nVYce3HH+Hcc91kpn6vhc6LHrhI\nwpk5aRLZ993Hc2vXHnruyZfWMrMhtO+mJF6eBg3cisyxY+HOO72OJjFpMyuRo5g8ZMhhyRvgubVr\nmfLaax5F5B933umGUSQ6lMBFjuLAjh1lPp9aUBDjSPync2fYssVVo0jkKYGLlCMYDDJ27Fi+37Wr\n7Ndr1oxxRP5zcDLz//7P60gSkxK4SDlSU1NJT0+n6/0v0iul6WGvPdG0KZ0HDvQoMn/p29eVE+bn\nex1J4lEVikgFrIUrr4SWjSdRe+NrpBYUEKxZk84DB2oCMwQ9ergdG/v18zoS/9B+4CKVtHv37jJr\nuEePhmHD4MsvIcnX5oTl00/hhRfcClapHCVwkaPIz88nJyeHdevWMWDAAKpVq3botU2b3D4n2dmQ\nnu5hkAmgqAgaNYJAAM45x+to/EEJXKQcwWCQefPmMWfOHC644ALat29PzSMmJW+6Cc480/UcJXyP\nPALVqsGf/+x1JP6gBC5Shh9++IEJEyZQr149MjMzy9zy9bPP4IEHYNkybYkaKcuXQ2Ym/Oc/kJrq\ndTTxTwlcpAy7du3ixx9/LHe3wF273NmOo0frbMdIa90ann3WJXKpmBK4SBX07+/+VO1y5A0fDjNn\nwj//6XUk8U8JXJKatZaCggKOCWEMZNIkGDDADZ1oA6bI274dGjeG9euhXj2vo4lvMTkTUyQe5eXl\nMWrUKKZPn17p92zd6vbuGD1ayTta6tVzwyfqgUeGeuCSUA6WBebm5tKxY0fS09MrfSrOzTe7HfRe\nfTXKQSa5rCz4wx9g/nyvI4lvGkKRpLJw4UKmTZtGenp6mWWBFfngA5dUlixR1Um0BYNw+ukwebI7\nFFrKpgQuSWXVqlXUr18/5JPgv/sOLrkEPvnEHcYr0ff4425xz0sveR1J/FICFzmKoiLIyIBf/cqd\nMi+xsXq1+7lv2KAtCsqjSUxJSPv37ydSnYKnn4batd0qQYmdc85xwyhTp3odib+FlcCNMTWNMfOM\nMV8ZY1YYY56PVGAiR7LWsnTpUoYOHcqGDRvCvt/UqfD22/Duu27faomtPn3gvfe8jsLfwh5CMcYc\na63da4xJA2YDD1trZ5d6XUMoEra8vDyysrIoLi6ma9euNGzYMKz7ff+9G/d+913o1ClCQUpINm+G\nZs3cMIrKNv9XTA41ttbuLXlYHUgFtoV7T5GDCgsLycrKqlJZYPn3hOuug7vvVvL20oknwuWXw8SJ\ncMstXkfjT5HogacAi4GmwAhr7e+OeF09cKmy4uJi5s6dy8UXX0yNGjXCvp+17lCBnTth3DgI83eB\nhOmDD9yWBVOmeB1J/IlVD7wYSDfG1AWyjTEZ1tpA6WsGDx586HFGRgYZGRnhNitJIiUlhXbt2kXs\nfsOHuwUkc+cqeceDX/0KfvtbyMuD007zOhpvBQIBAoFASO+JaBmhMeb3wD5r7culnlMPXCqlqKiI\ntCjWlGVnu4/qX3wBTZse/XqJjX79oHlzePhhryOJL1EvIzTG1DfGHF/y+BigM7AknHtK8iksLGTq\n1Km8+eabFBcXR6WNRYugd2/48EMl73jTu7eqUaoq3O7OKcA7JePgKcB71tqc8MOSZGCtZdmyZeTk\n5NCkSRP69OlDShTq+datcx/V33wTIjgaIxHSvj3s2OF2gGzZ0uto/EUrMcUTP/74I5MmTYpYWWB5\nNm6EK66A+++He++NShMSAU88AQcOaGl9aVpKL3Frw4YNbNmyJSJlgeXJy3PLte+8Uyst492KFdC5\ns9uXRsetOUrgkrTy8lzPu18/+N3vjn69eK9VK3jxRbjySq8jiQ/aC0XiQrQmJsuzcqUb6+7fX8nb\nT7S0PnTqgUvUbNu2jezsbBo0aMAVV1wRkzbnzIGePeEvf4Fbb41JkxIhP/7oygnz8qBWLa+j8Z56\n4FKhaHWMD5YFjhw5ktNPP53LL788Og2VYq2rMunRw+1vouTtPw0awKWXwscfex2Jf6gHnsAOHIAF\nC9zXokWQm+t6Nz/9BPv3uwReowbUrQs/+xk0aeJqpM8+241HpqeHfjrNsmXLmDp1Kk2aNKFTp07U\nicEuRXv3ugqTBQtgwgS3Van40z/+4YZRPv/c60i8p0nMJFRUBP/6lzs0NisLzjgD2rZ1Cfmcc6Bh\nQzjpJKhZ022hWlDg9gXZssXVS69dC6tWwcKFbiz57LOhQwe36VOHDi7ZA8ycNInJQ4aQVlhIUY0a\nXDVoEO27dWPWrFk0btw4amWBR5o9G/r2dd/jiBFub2/xr7173ZL6lStdjzyZKYEnkZ07YehQeOMN\n9x/g1lvd4pVw9pcoKIClS2H6dMjJgS+/hPPPh/PPmESd2ffx17y1h659smlTMv/+d9p36xaB7+bo\ntmyBwYPdysphw9zQiSSGW291n/4eeMDrSLxVmQSOtTaqX64JiZa9e619/nlr69e3tk8fa5csiV5b\n+/ZZO3WqtTeceZW1btj5sK+nMjOj13iJ3butffFF9/3ee6+1W7ZEvUmJsSlTrL3wQq+j8F5J7qww\nv2oS08cmT3Y94gULYOZMN3mXnh699mrWdEMpjU/aUebrC2cUMGAAfPop5OdHtu31692ZlWee6b7f\n2bPdJ44TTohsO+K9K65w8zTLl3sdSfzTcaI+tHs3DBjgkvbw4dC1a2zaPVgW+P2uXWW+fvaFNTm1\nEbz6Ktx8M1x8sVuUcXBC9JRTKt/WgQOwZAkEAjB+PHz7LfzmNy55N24ckW9H4lRqKvTq5SYzX3jB\n62jim8bAfWbxYrjxRrdE/G9/i0297P79+5k5cyaLFy+mXbt27N+6lakPPshza/87Bv5E06Z0KTUG\nvmePS76BgEvES5a408ebNnWH2TZsCMcd56pc0tJcj33XLne8Vm6uO7W8cWN3YkuPHu771enlyePr\nr+GXv4T//Cd5zyvVJGaCef99N7Hz2mtw002xa3fPnj1MmzaNjIyMQ2WBMydNYsprr5FaUECwZk06\nDxxY4QSmte4cyvXr3X4XGza4pF1Q4HrbtWu7cxEbNoSzznIVM8cfH6vvUOJRerr7NBejNWBxRwk8\nQVgLTz/tPlJ+9hmcd57XEYlE3yuvuHHwUaO8jsQbSuAJIBh0GzKtWgUffQQnnxzd9qy1UdsdUCQU\nP/wALVq4xWfHHut1NLGnpfQ+d+CAm7jbuNHVYUczeQeDQb744gvGjh0bvUZEQnDKKdC6NXzyideR\nxC8l8Di1fz/8+tdunPjTT6PbA8nNzWXEiBGsX7+eq666KnoNiYRIOxRWTEMocSgYdD3vfftg3Dio\nXj067RwsC9yyZQuZmZk0a9ZMwycSV/bscauJV6+O/vBhvNEYuA9Z62q8V6xwG/rUrBm9tr766iv2\n7NlDmzZtonoavEg4brnFrSW47z6vI4mtqCdwY0wj4F3gJMACb1prhxxxjRJ4CAYPdkMm06e7OmmR\nZDdlCjz+uNtgLZnEYhLzAPCAtfY8oC1wrzHm3DDvmbTee88th//8cyVvkYM6dnQVKStXeh1J/Akr\ngVtrf7TWflXyOB9YCZwaicCSzZdfwoMPut73SSdF7r75+fl8/PHHLFmyJHI3FYmh0kvr5XARq0Ix\nxpwJXAjMi9Q9k8WGDXDddTB6dOQW6RwsCxwxYgTHHHMMLVq0iMyNRTzQpw+MGRO9U6T8KiIzV8aY\n2sB44L6SnvhhBg8efOhxRkYGGRkZkWg2IezbB9dc45bId+8emXvm5uaSlZVFvXr16Nu3L/Xr14/M\njUU80rKl21ph5ky3L04iCgQCBAKBkN4TdhWKMaYa8BnwubX2b2W8rknMCtx1F+zY4U7QiUQFn7WW\nCRMm0LJlS5UFSkJ5+WU3Dv7WW15HEhuxqEIxwDvAVmttmednKIGX7x//cFUnCxdq0lLkaL7/3u1/\nn5cX+lmtfhSLKpR2QG/gCmPMkpKvLmHeMymsXu3qWj/4QMlbpDJOPdXtMa+l9f8VbhXKbGttirU2\n3Vp7YclXVqSCS1T79sENN8Bzz1X9BJ28vDzGjBnD3r17IxucSBzT0vrDaSWmBwYNckdGjR0b+rh3\nfn4+OTk55Obm0rFjR9LT0zXOLUkjP9/tGf/vf0e23DYeaSl9HJoyBe64w532Xq9e5d8XDAaZN28e\nc+bM4YILLqB9+/bUjOY6e5E41acPXHKJ6wglMiXwOLN9uyuHGj3anRUZih9++IFp06aRmZmpskBJ\napMnw5NPuvNRE5kSeJzp1Qvq14chQ45+rYiULRh056pOnQrnJvDGHTrQIY7885/uQGKdsi0SntRU\n6N0b3n7b60i8px54DGza5IZOPvvMjd2Vx1rLsmXL2Lx5M1eGOsYikkRWr3YrMr/7DqpV8zqa6FAP\nPE4MHAh9+1acvPPy8hg1ahTz58+nefPmsQtOxIfOOQeaNnU7dyYz7eIfZR9/DF99Be+8U/brKgsU\nqZrbb3cn1l99tdeReEdDKFG0Y4db+jtmDHToUPY1kydPBlBZoEiIdu92k5krV0KDBl5HE3mqQvHY\nnXdCSgq8/nr511hr1eMWqaI77oDmzeGRR7yOJPKUwD0UCLgFB998A3Xreh2NSGKaMwf69XNnyCZa\nP0iTmB7Ztw/694fhw13yLiwsZOrUqWzYsMHr0EQSymWXuYPAv/zS60i8oQQeBS+8ABdcAN27W5Yu\nXcqwYcPYvXs3xx9/vNehiSQUY9xkZrLsEX4kDaFEyMxJk5g8ZAj7txcye0kN+r/Si+JaRRQXF9O1\na1caNmzodYgiCemHH6BFC3c0Ye3aXkcTORoDj5GZkyaRfd99PLd27aHnBpx4Iq2eeorbBg7UJKVI\nlF19NfTsCbfd5nUkkaMx8BiZPGTIYckbYOjmzaz717+UvEVi4PbbYeRIr6OIPSXwCEgrLCzz+dSC\nghhHIpKcuneH9evh66+9jiS2lMDDsHXrVj7//HMO1KhR5utBLcwRiYm0NFf5VdGai0SkBF4FhYWF\nTJkyhbfeeos6derQrPu99Eppetg1TzRtSueBAz2KUCT59OvnTrnKz/c6ktgJexLTGDMK6Ab8ZK39\neRmvJ8wk5sHdAnNycmjSpAmdOnWiVq06tG8P7VpOovq610gtKCBYsyadBw6kfbduXocsklR69ICu\nXd0qaL+LSRWKMeZyIB94N9ET+KpVq5g1a9ZhZYFvv+0W7Myd6/YpFhHvZGfDY4+5vff9Xj8QszJC\nY8yZwKeJnsAPfh8HK0u2bXP1p5MmQatWXkYmIgDFxXD22fD++9C2rdfRhEdlhBFmjDmsLPCJJ+D6\n65W8ReJFSgrcdVfyTGbGZD/wwYMHH3qckZFBRkZGLJqtstzcXPbu3UvLli3LvWb+fLfX98qVMQxM\nRI7qttugWTP3CflnP/M6msoLBAIEAoGQ3qMhlFK2bdtGdnY2W7ZsoWvXrpx11lllXhcMQuvWcP/9\nbsdBEYkvffq4/YgeftjrSKpOY+CVVFhYyKxZs1i8eDHt2rWjTZs2pKWV/+Fk2DAYNw6mT/f/RIlI\nIlq4EK67DtaudTXifhSrKpSxQAfgBOAn4A/W2tGlXo/7BD5u3DiqVatGp06dqFOnToXXbtrkTtkJ\nBOC882ITn4iE7he/cJ+Sr7/e60iqRptZVVJRUVGFPe7Sbr0VTj4Z/vKXKAclImGZMAH++ld36IMf\nKYFH2MyZ0Lu3O/0jkbatFElERUVw1lnwwQduzspvVEZYSjAYZO7cuezevbtK7z9wAO65x/1GV/IW\niX9paTBoEPztb15HEj1JkcDXrFnDiBEjWLduHcFgsEr3GDIETjvNTYyIiD/ccQdkZcHGjV5HEh0J\nPYRSuiwwMzOTs88+u0r3yctzJUlz57r6UhHxj/vug2OPheef9zqS0CT1GPjevXsZPnw4l1566VHL\nAo/mxhvd8txnn41ggCISE+vXwyWXuJLCunW9jqbykjqBA+zfv5/q1auHdY8pU9zS3OXL4ZhjIhSY\niMRUnz6u7Pexx7yOpPKSPoGHq7AQWraEV15xJ36IiD8tXw6dOsG6dW44xQ+SogolPz+fRYsWReXe\nzz8P556r5C3id+edB5ddBm+95XUkkeXbHngwGGTevHnMmTOHCy64gM6dO0f0AOEVK6BDB/jqK1d9\nIiL+tmCBqyLLzYUwR1ZjImGHUHJzc8nKyqJevXpkZmZSv379iN6/uBguvxx+8xtX+y0iieGqq+Cm\nm9wp9vEuIRP44sWLmTNnDpmZmTRr1iyive6DRoxwG8LPmuX2FxaRxDBzJvTtC6tWQbVqXkdTsYRM\n4Pv37yclJSWsssCK5OVBejrMmOFO2xGRxJKZCT17uuqyeJaQCTzaevRwlSfPPON1JCISDQsXwrXX\nwpo18V0a7OsqlLy8PDbGeP3rhAnuo9UTT8S0WRGJoYsvhjZt3L7+fhd3PfD8/HxycnLIzc2lW7du\nNG/ePIrR/demTW65/MSJcOmlMWlSRDyyYgVkZLiKlOOO8zqasvlqCOXIssAOHTpQo0aNqMZ2kLVu\nTKx5c//tlyAiVXPbba5E+LnnvI6kbL5K4O+++y6pqalRKQs8etvw8suuTjRGvzNExGMHN6lbsAAa\nN/Y6mv/VdjWwAAAH2klEQVTlqwSen59PrVq1olIWWJENG6BVK5g82VWfiEjyeO45WLIExo/3OpL/\nFZNJTGNMF2PMKmPMGmPMo1W9T+3atWOevINB9zFq0CAlb5Fk9OCDsGiRO6Dcj8JK4MaYVGAo0AVo\nAdxsjDm3vOuttSxfvrzKhypE2vPPuyT++ONeRyIiXjjmGDd8OnAg7N/vdTShC3c1TGsg11r7LYAx\n5p/ANcDK0hc9lZnJRb16sb2oiOLiYho1asRxHk/9zpwJQ4e6376pqZ6GIiIe6tnTbXL10kvw5JNe\nRxOacIdQTgM2lPr7xpLnDvOnyZOZ+tBDpOzZQ79+/TxP3lu2uH1ORo3SRlUiyc4Yt33Gq6/C6tVe\nRxOacBN4pWdAh2/dyrp//Svm49xHCgbd5u433QS//KWnoYhInDjjDPj97+HOO91mdn4R7hBKHtCo\n1N8b4Xrhhxlc8uesVasIBAJkZGSE2WzVPf64O6jhz3/2LAQRiUMDBsDPfubWhXghEAgQCARCek9Y\nZYTGmDRgNdAJ+B6YD9xsrV1Z6ppDLfw+M5Nns7Kq3F643n8fnn4a5s+HE07wLAwRkaOqTBlhWD1w\na22RMWYAkA2kAm+VTt6lDTixKb8eODCc5sIybx488IArF1LyFpFEEPaerNbaz4HPK7rmgXaZTFw2\nkNtO7hZuc1WyahVccw2MHg3nn+9JCCIiERezlZgTJ7oFM/PnwymnRLXJw2zYAL/4Bfzxj3DrrbFr\nV0QkHHG1nWyPHm4D9Wuvhfz82LT5/ffQubMr0lfyFpFEE9O9UKx1ZTpr18KkSdHdTP0//4FOnaBf\nP3jssei1IyISDXHVAwdXMP/6624IpWdP2LcvOu0sWwbt27uyICVvEUlUMT+RJzUV3nnHVYJ06uRW\nRUbSJ5+4+77wAtx/f2TvLSISTzw5Ui0tze3B3aEDtG3r9iMJV0GB21nsnnvgs8/g5pvDv6eISDzz\n7EzMlBS3G+Bzz0HXru5xYWHV7jVzJlxyCXz3HSxd6s67ExFJdJ4fanzjjW6RzZw5rkZ77Fg4cODo\n77PWvefaa+GWW+Cpp2DcOC3SEZHkETcn8gBkZ8OLL7odwXr2dCWALVrAqae6hL19O6xcCTNmwMcf\nu2GTgQOhf//oVrSIiMSar45UK23FCvj0U5g2DdascfXcKSlQt647eLhNG+jWDdq1c8+LiCQa3yZw\nEZFkF3d14CIiEjlK4CIiPqUELiLiU0rgIiI+pQQuIuJTSuAiIj6lBC4i4lNK4CIiPlXlBG6MucEY\ns9wYEzTGXBTJoERE5OjC6YF/DfQAZkYolrgUCAS8DiEsfo7fz7GD4vea3+OvjConcGvtKmvtvyMZ\nTDzy+z8CP8fv59hB8XvN7/FXhsbARUR8Kq2iF40xU4AGZbz0hLX20+iEJCIilRH2boTGmOnAQ9ba\nxeW8rq0IRUSq4Gi7EVbYAw9BuY0cLQAREamacMoIexhjNgBtgUnGmM8jF5aIiBxN1A90EBGR6Iha\nFYoxposxZpUxZo0x5tFotRMtxphRxphNxpivvY4lVMaYRsaY6SULrb4xxgzyOqZQGGNqGmPmGWO+\nMsasMMY873VMVWGMSTXGLDHG+G7C3xjzrTFmWUn8872OJxTGmOONMeONMStL/v209TqmyjLGnFPy\nMz/4tbOi/79R6YEbY1KB1cCVQB6wALjZWrsy4o1FiTHmciAfeNda+3Ov4wmFMaYB0MBa+5Uxpjaw\nCLjWZz//Y621e40xacBs4GFr7Wyv4wqFMeZBoBVQx1p7tdfxhMIYsx5oZa3d5nUsoTLGvAPMsNaO\nKvn3U8tau9PruEJljEnB5c/W1toNZV0TrR54ayDXWvuttfYA8E/gmii1FRXW2lnAdq/jqApr7Y/W\n2q9KHucDK4FTvY0qNNbavSUPqwOpgK8SiTGmIfBLYCQVTPLHOd/FbYypC1xurR0FYK0t8mPyLnEl\nsLa85A3RS+CnAaUb3VjynMSYMeZM4EJgnreRhMYYk2KM+QrYBEy31q7wOqYQvQo8AhR7HUgVWWCq\nMWahMaa/18GEoDGw2Rgz2hiz2Bjzf8aYY70OqopuAv5R0QXRSuCaGY0DJcMn44H7SnrivmGtLbbW\npgMNgfbGmAyPQ6o0Y0x34Cdr7RJ82Ist0c5aeyHQFbi3ZEjRD9KAi4Dh1tqLgD3AY96GFDpjTHXg\nV8C4iq6LVgLPAxqV+nsjXC9cYsQYUw2YALxvrf3I63iqquTj7yTgYq9jCcFlwNUl48hjgY7GmHc9\njikk1tofSv7cDEzEDYv6wUZgo7V2Qcnfx+MSut90BRaV/PzLFa0EvhBoZow5s+Q3yY3AJ1FqS45g\njDHAW8AKa+3fvI4nVMaY+saY40seHwN0BpZ4G1XlWWufsNY2stY2xn0MnmatvcXruCrLGHOsMaZO\nyeNawFW43UfjnrX2R2CDMebskqeuBJZ7GFJV3Yz75V+hSK3EPIy1tsgYMwDIxk1AveWnCggAY8xY\noANwQsmCpT9Ya0d7HFZltQN6A8uMMQcT3+PW2iwPYwrFKcA7JbPwKcB71tocj2MKh9+GFE8GJrp+\nAGnAGGvtZG9DCslAYExJ53Et0NfjeEJS8kvzSuCocw9ayCMi4lPaTlZExKeUwEVEfEoJXETEp5TA\nRUR8SglcRMSnlMBFRHxKCVxExKeUwEVEfOr/AT5at47r+ZddAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "metadata": {}, "input": ["pt.plot(plot_x, interpolant)\n", "pt.plot(plot_x, f(plot_x), \"--\", color=\"gray\")\n", "pt.plot(x, f(x), \"or\")"], "prompt_number": 58, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["* For $f(x)=x$, why does it do what it does at the interval end?\n", "* What happens if we increase the number of points?\n", "* What if the function is actually periodic (e.g. the function with `abs` above?)?\n", "* What if the function has jumps?"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["# Answers\n", "#\n", "# * Because we're interpolating with periodic functions--so the interpolant is forced to be periodic.\n", "# * We observe a distinct \"overshoot\". This overshoot is called \"Gibbs phenomenon\".\n", "# * Periodic: no Gibbs at interval end.\n", "# * Gibbs can also happen in the middle of the interval."], "prompt_number": 59, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["## Computing the coefficients"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["B = V.T.dot(V)"], "prompt_number": 60, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [], "metadata": {}, "input": ["B[np.abs(B)<1e-12] = 0"], "prompt_number": 61, "collapsed": false, "language": "python", "cell_type": "code"}, {"outputs": [{"prompt_number": 62, "text": ["array([[ 5. ,  0. ,  0. ,  0. ,  0. ],\n", "       [ 0. ,  2.5,  0. ,  0. ,  0. ],\n", "       [ 0. ,  0. ,  2.5,  0. ,  0. ],\n", "       [ 0. ,  0. ,  0. ,  2.5,  0. ],\n", "       [ 0. ,  0. ,  0. ,  0. ,  2.5]])"], "metadata": {}, "output_type": "pyout"}], "metadata": {}, "input": ["B"], "prompt_number": 62, "collapsed": false, "language": "python", "cell_type": "code"}, {"source": ["* What do you observe?\n", "* How could this be useful for computing the inverse?\n", "* What is the normalization?\n", "* This is a pretty special matrix. What is the cost of computing $Ax$ with it?\n", "* The so-called [Fast Fourier transform](https://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT) computes a version of $Ax$ (with complex exponentials) in $O(n\\log n)$!"], "metadata": {}, "cell_type": "markdown"}, {"outputs": [], "metadata": {}, "input": ["# Answers:\n", "# \n", "# * V's columns are orthogonal. (though not normalized)\n", "# * The transpose of $V$ (with appropriate normalization) its inverse. This makes Fourier coefficients cheap to compute.\n", "# * The normalization is n for the first entry, and n/2 for all the ones after that.\n", "# * Computing Ax costs n**2 operations."], "prompt_number": 19, "collapsed": false, "language": "python", "cell_type": "code"}]}], "nbformat": 3}