{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\bkt}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\dsum}[1]{\\displaystyle\\sum}$\n", "$\\newcommand{\\spade}{\\bkt{\\spadesuit}}$\n", "$\\newcommand{\\club}{\\bkt{\\clubsuit}}$\n", "\n", "Polynomial Interpolation\n", "==\n", "\n", "1.1 **Introduction**\n", "\n", "Given the values of a function $f(x)$ at $n+1$ distinct locations of $x$, say $\\{x_i\\}_{i=0}^n$, we could approximate $f$ by a polynomial function $p_n(x)$ of degree $n$ that satisfies\n", "\n", "$$p_n\\bkt{x_i} = f\\bkt{x_i}$$\n", "\n", "We can construct the polynomial $p_n(x)$ as $p_n(x) = a_0 + a_1 x + a_2 x^2 + \\cdots + a_n x^n$. The $n+1$ coefficients are determined by forcing $p_n(x)$ to pass through the data points. This leads to $n+1$ equations in $n+1$ unknowns, $a_0,a_1,\\ldots,a_n$, i.e.,\n", "$$y_i = a_0 + a_1 x_i + a_2 x_i^2 + \\cdots + a_n x_i^n$$\n", "for $i \\in \\{0,1,2,\\ldots,n\\}$. This procedure for finding the coefficients of the polynomial is not very attractive. It involves solving a linear system, whose matrix is extremely ill-conditioned. See below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import scipy as sp;\n", "import numpy as np;\n", "\n", "from scipy.stats import binom;\n", "import matplotlib.pylab as pl;\n", "from scipy import linalg\n", "from numpy import linalg\n", "from ipywidgets import interact;" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdXZ/vHvQ5SiglgRrCVKUEGJzAQERRSkEvpSJm1F0ReEglZEwIGhip7DYEHqT3kVBEJIBC3RFkUQZFRBKSoJgiRBa6RYAgoUGZTKmPX7IyENIUDgnGSf4f5cVy85a+9kP9duvbt41j5rm3MOERGJHhW8LkBERMqXgl9EJMoo+EVEooyCX0Qkyij4RUSijIJfRCTKKPhFRKKMgl9EJMoo+EVEosw5XhdQkksuucTFxcV5XYaISFjJyMj4t3Ou+unOC8ngj4uLIz093esyRETCipl9U5rz1OoREYkyCn4RkSij4BcRiTIKfhGRKKPgFxGJMgp+EZEQ4vP5yvwaCn4RkRDi9/vL/BoKfhGRELF06dJyuY6CX0TEYz6fDzPjtttuA8DMMLMya/uE5Dd3RUSihXOOI0eOANChQwcWL16Mc65Mr6kZv4iIRw4ePMg999zD2LFj6devH/Pnzy+X62rGLyLigV27dtGtWzc+/PBDxo0bx9ChQzEznn766TK/tmb8IiLlyOfz8fXXX3PDDTfw6aefkpaWxrBhwzCzwuNlzcq6l3Q2EhISnHbnFJFIZGZccskl5OXl8fbbb9O6detg/u4M51zC6c5Tq0dEpBxs3ryZv/3tbwBcdNFFLFy4kDp16nhSi1o9IiJBdKxV45wjIyODp556iksvvZTatWvz+OOPA5CTk0PdunXLpa1TkqDP+M3sSuAJoKpz7o6CsZuAngXXi3fO3RDs64qIeO3w4cP4/X527tzJvHnzyM3NpUKFCtx4440MHTqUzp07U7du3TJ/XPN0SjXjN7MZZrbDzDKLjSea2ZdmlmNmwwGcc5ucc32Lnuec+9A59wDwDvBKsIoXEQkVWVlZNG7cGIDU1FQSEhJISUnhu+++Y+XKlTz66KOetXaKK22rJxVILDpgZjHAJKAjEA/cZWbxp/k9dwOzz7BGEZGQ5Zyjc+fO1K9fn+zsbAD+85//MHfuXDZv3kz16se/Arc8Htc8nVK1epxzK80srthwCyDHObcJwMzSgC5Adkm/w8yuAPY65/addbUiIiHkhx9+4IEHHmD+/PnceuutvPrqq1x22WWnbOV41dcvKpDF3ZrAliKfc4GaZlbNzKYATcxsRJHjfYGUk/0yM+tvZulmlr5z584AyhIRKXvr1q2jWbNmpKWlMXr0aBYvXswvfvELr8sqlUAWd62EMeec2wU8UMKBU/79xjk3DZgG+c/xB1CXiEiZefrpp7n00kt55JFHqFatGu+//z5t2rQ57nioCyT4c4HLi3yOBbYFVo6ISOjas2cPo0aNAqBjx4688sorJ/TwQ6GVczqBtHrWAHXMrLaZVQR6APOCU5aISGhZvnw5DRs2BODZZ5/lnXfeOSH0w0VpH+ecDawGrjGzXDPr65w7AjwELAY2Am8457LKrlQRkfK3f/9+WrRoQfv27dmyJX9Zc+jQocTExITF7L4kpX2q566TjC8EFga1IhGRELFq1Sp69erF119/zeDBgxk7diwXXHCB51/ACpS2bBARKcLn83HgwAGGDh3KTTfdRF5eHh988AHPP/88559/vtflBYU2aRMRKcLv9/PXv/6V7Oxs7r//fiZMmECVKlUKj4fDUzuno+AXESF/n52xY8cCsHfvXhYtWkSHDh1OOC9c+/pFqdUjIlFvwIABVKxYEb/fD8DWrVtJTEyMiJAviWb8IhK18vLymDhxIsnJyVxyySVMnTqV22+/PewXb09HwS8iUWnz5s307t2bFStW0LlzZ6ZNm8all17qdVnlQsEvIlHl6aefplatWgwePBiAlJQUevXqVfjO20hYvD0dvXNXRKLG9u3bCzdSu+WWW0hNTaVWrVoeVxU8pX3nrhZ3RSTiOeeYPXs21113HQDPP/88y5cvj6jQPxMKfhGJaN999x316tXj7rvvZteuXQAMGTIkrLdcCJSCX0QiknOO1157jfj4eDZv3syECRM4cuRI4THnnIJfRCQS+Hw+vv32W7p27co999zDtddey7p163jssceIiYnxuryQoOAXkYjhnMPv9xMfH8+SJUt47rnn+PDDD7n22msLz4mGp3ZOR49zikhE2Lp1Kw88kP/yv+uuu44ZM2ZQt27dE86L1vZOUZrxi0hYc87RuXNnYmNjeeedd4D87ZSvueYahfxJaMYvImHrm2++oX///ixZsoSbb76Z6dOnU6dOnYjfciFQmvGLSFjx+Xzk5eXx8ssvU79+fVatWsWkSZN47733uPrqq70uLyxoxi8iYcXv97NixQo++OAD2rdvT1JSEnFxcYXHtXh7egp+EQkLeXl5vPjiiwCsXbuWpKQk+vbtW7jHzjHq65+eWj0iEvIGDhxITExM4cZq+/bto1+/foX758uZ0YxfRELW0aNHefHFF0lOTqZq1apMnDiR3r17a/E2QJrxi0hIOdaq+cc//sHNN9/MkCFDaNeuHVlZWfTq1cvb4iKEZvwiElL8fj8XXnghTzzxBJUqVWLmzJncc889UbVfflkL+n78ZnYl8ARQ1Tl3R5HxC4CVwNPOuXdO9Tu0H79IdPryyy8Lt1fo1KkTU6dO5Ze//KXHVYWPoO7Hb2YzzGyHmWUWG080sy/NLMfMhgM45zY55/qW8GuGAW+U5noiEl2eeuopzOy4PXXeeecdpk2b5mFVkau0rZ5U4CVg5rEBM4sBJgG/AnKBNWY2zzmXXfyHzaw9kA1UCrRgEYksWVlZLFq0CIAuXbrw9ttva/G2jJVqxu+cWwl8X2y4BZBTMMM/BKQBXU7yK9oCLYG7gX5mpkVlkSh3+PBhxowZQ5MmTfjnP/9JWloab731ltdlRYVAArgmsKXI51ygpplVM7MpQBMzGwHgnHvCOTcY+AuQ5JzLK/7LzKy/maWbWfrOnTsDKEtEQpnP52PdunW0aNGCkSNH0r17d7Kzs7nzzjsxMy3eloNSL+6aWRzwjnOufsHn3wIdnHO/L/h8L9DCOTcw0KK0uCsSmQ4ePEilSpU455xzqFatGi+//DLdunXzuqyIUR4vW88FLi/yORbYFsDvE5EIlp6eTrNmzQC4++67yc7OVuh7JJDgXwPUMbPaZlYR6AHMC05ZIhIpDh48yE033UTz5s3JysoCYObMmVSrVk376niktI9zzgZWA9eYWa6Z9XXOHQEeAhYDG4E3nHNZZVeqiISbTz/9lKZNm/LRRx/Rp08fdu/eDehl514r7VM9dznnLnPOneuci3XOJReML3TO1XXOXeWcG1u2pYpIOPD5fBw4cIDhw4fTqlUr9u3bx7vvvktycjIXXXSR1+UJ2rJBRILM7/fz+uuv88UXX/D73/+eP//5z1StWrXwuJ7a8Z6CX0SC4qeffioM9f3797No0SI6dOhwwnlq73hPX6QSkYD16dOH888/nwkTJgCwZcsWEhMTFfIhSjN+ETlr+/fvZ8SIEaSmphIXF8f06dNp3769tlwIcZrxi8gZOTaLf++992jQoAEvvfQSAwcOZMOGDdx6663eFielohm/iJwRv9/Pt99+y7Rp06hTpw4rV66kdevWhce1eBv6NOMXkVI7tovm9OnTefzxx1m/fv1xoQ9avA0HCn4ROa3hw4djZnTs2BGAvLw8JkyYwPjx4z2uTM6GWj0ickoLFixg1qxZxMTEMHz4cMaOHavF2zCnGb+IlGj37t307t2bTp06cfHFF/PJJ58wZswYr8uSINCMX0SO4/P5aNasGffffz87duxg5MiRPPnkk1SsWBHQ4m0kCPrL1oNB+/GLeOP777+nWrVqADRs2JCUlBSaNm3qcVVSWuWxH7+IRJC33nqL+Ph4IH9Wv2bNGoV+hFLwi0S5HTt2cN1119G9e3e2b98O5D+r/7Of/UyPZkYoBb9IlHLOMXv2bOLj48nJyWHs2LEcOnSo8Jj2y49cCn6RKOPz+di2bRtdu3bl7rvv5uqrr+azzz7jj3/8I+eee67X5Uk5UPCLRBHnHH6/n/j4eJYsWcJzzz3HqlWrCnv7oKd2ooGCXyRK/Otf/yIxMRGARo0a8fnnn/PII48QExNz3Hlq70Q+Bb9IhHPO8Zvf/IZatWqxZMkSAFauXEndunUV8lFKX+ASiWCbN2+mX79+LFu2jHbt2pGcnEzt2rW15UKU04xfJML4fD7y8vJ4+eWXadCgAR9//DFTp05l2bJlxMXFeV2ehADN+EUijN/vZ8WKFXzwwQf86le/IikpiVq1ahUe1+KtKPhFIkReXh6TJ08GICMjg6SkJPr27YuZHXee+voS9FaPmV1pZslm9rdTjYlI8AwcOJCYmBgGDhwIwA8//EC/fv3w+/0eVyahqFTBb2YzzGyHmWUWG080sy/NLMfMhgM45zY55/oWPa+kMREJ3JEjR5gwYQLTp0/noosuIjU1FdA3b+XUSjvjTwUSiw6YWQwwCegIxAN3mVn8iT8qIsF0LMw3bNhAq1atGDp0KB06dCA7O5tevXp5W5yEhVL1+J1zK80srthwCyDHObcJwMzSgC5AdjALFJHj+f1+KlSowJgxY6hatSppaWn87ne/K+zla/FWTieQHn9NYEuRz7lATTOrZmZTgCZmNgKgpLHizKy/maWbWfrOnTsDKEskcmVkZAD54X7HHXeQnZ3NnXfeedwCrto7cjqBBL+VMOacc7uccw84565yzv2pYPCEsRJ+cJpzLsE5l1C9evUAyhKJPE8++SRmRkLCf9+xMXv2bCZNmuRhVRKuAnmcMxe4vMjnWGBbYOWISHGffvopb731FgC9e/cmNTVV37yVgAQy418D1DGz2mZWEegBzAtOWSJy4MABhg0bRqtWrdi7dy8LFy4kJSXF67IkApT2cc7ZwGrgGjPLNbO+zrkjwEPAYmAj8IZzLqvsShWJDj6fj9WrV9OkSROeffZZ+vTpQ1ZWFh07dgS0eCuB08vWRULITz/9xPnnn4+ZERsby/Tp07ntttu8LkvChF62LhJmPvzwQxo1agRA//79yczMVOhLmVDwi3jshx9+oEWLFrRp04avvvoKgKlTp1K1alU9millQsEv4qElS5ZQv3590tPTGTRoED/++COgLRekbCn4RcqZz+dj9+7d9OnThw4dOnDeeefx0Ucf8cILL3DBBRd4XZ5EAQW/SDk79rLzmTNnMmLECNatW8cNN9xQeFxP7UhZ0378IuVk586dhdsm16hRgwULFtC0adMTzlN7R8qaZvwi5eC3v/0tNWrU4PXXXwfg888/p1mzZgp58YRm/CJlaPv27QwYMIA5c+bQrFkzUlJSaNiwobZcEE9pxi8SZD6fD+ccaWlpXHfddcyfP58//elPfPzxxzRo0MDr8kQ04xcJNr/fz/r165k7dy4tWrQgJSWF+Pj/vqNIi7fiNc34RYLEOcerr74KwLvvvsuzzz7LqlWrjgt90OKteE/BLxIEQ4YMoUKFCtx7770AHDx4kKFDhzJmzBiPKxM5kVo9IgHIy8sjKSmJ5ORkzj//fJ555hkGDx6sxVsJaZrxi5yhY62anJwcbr31Vh544AFatGjBhg0bGDRokLfFiZSCZvwiZ8jv91OlShVGjhzJueeeS1JSEn379tXLziVsaD9+kTOQlZVF/fr1AejcuTOTJ0+mZs2aHlclkk/78YsE0ciRIzGzwtAHmDdvHklJSR5WJXJ21OoROY3169ezYMECAO68805ef/11Ld5KWNOMX+QkDh06hM/nIyEhgW3btvHmm2+SlpbmdVkiAVPwixTj8/lYu3YtzZs3x+/306NHD7KysujWrRugxVsJf1rcFSni4MGDVKpUiZiYGGrUqMHUqVP5zW9+43VZIqWixV2RM/TJJ5/QrFkzAO69916ysrIU+hKRFPwS9fbv30+rVq1o2bIlWVlZAKSmpnLxxRdrXx2JSOUS/GZ2hZnNM7MZZja8PK4pUhrLly+nQYMGfPzxx/zhD39g7969gF52LpHtrIO/IMR3mFlmsfFEM/vSzHKKhHxdYIFzrg8Qf8IvEylHPp+PPXv20K9fP9q3b88555zDihUrmDx5MhdeeKHX5YmUuUBm/KlAYtEBM4sBJgEdyQ/4u8wsHvgM6GFm7wHvB3BNkYAde9n5jBkzGDp0KOvXr6dNmzaFx/XUjkS6s/4Cl3NupZnFFRtuAeQ45zYBmFka0AU4DDxd8DN/A1LO9roiZ6voy84vueQS5s2bR0LCiQ9AqL0jkS7YPf6awJYin3MLxhYBD5vZFGBzST9oZv3NLN3M0nfu3BnksiSaOedOeNn5hg0baN68uUJeolKwt2ywEsaccy4TuONUP+icmwZMg/zn+INcl0Sp7du38+CDD/Lmm2+SkJBASkoKDRo00JYLEtWCPePPBS4v8jkW2Bbka4ic0rGXnb/22mvEx8ezYMECxo0bx+rVq4/bZE0kWgV7xr8GqGNmtYGtQA/g7iBfQ+SU/H4/a9euZf78+bRs2ZIZM2ZQr169wuNavJVoF8jjnLOB1cA1ZpZrZn2dc0eAh4DFwEbgDedcVnBKFTk15xypqakALF26lOeee46PPvrouNAHLd6KBPJUz10nGV8ILDzrikTOwuDBg5k4cWLh5wMHDvDoo4+yb98+Bb1IMdqPX8La0aNHmTRpEtOnT6dy5cqMHz+eAQMGaPFW5BS0V4+EnWMz+OzsbFq3bs2gQYNo06YNWVlZPPjgg94WJxIGNOOXsOP3+4mJiWHMmDFUqVKFWbNm0bNnT73sXKSUtB+/hJU1a9bQokULAHr06MHEiROpUaOGx1WJhAbtxy8R5YknnsDMCkMfIC0tjcmTJ3tYlUh4UqtHQt7f//535syZA0C/fv1ISkrS4q1IADTjl5D1n//8hyFDhtC6dWsOHDjA0qVLmTZtmtdliYQ9zfgl5Ph8Ptq2bUvfvn35+uuvefDBBxk3bhxVqlQBtHgrEigt7kpI+fHHHwsD/sorryQ5OZlbbrnF26JEwoQWdyXsLF68uHATtUGDBvH5558r9EXKgIJfPLdr1y4aNWpEYmIi33zzDQATJ06kcuXK2m5BpAwo+MUzzjlef/116tWrR3Z2Nk8++SQ//fRT4TG97FykbCj4pdz5fD5yc3Pp0qULPXr0IC4ujoyMDEaPHk2lSpW8Lk8k4in4pVzl5eUVvux82bJlPPfcc6xevZqGDRsWnqOndkTKloJfys1XX31F27ZtAbj++uvJzMzkkUceISYm5rjz1N4RKVsKfilzR48e5bbbbqNu3bqsXLkSgGXLlnHVVVcp5EU8oOCXMpWdnc2NN97I0qVL6dy5M1u3bgW0eCviJQW/BJ3P5+Pw4cM888wzNGnShJycHP7yl78wd+5cfvnLX3pdnkjU05YNEnR+v5/58+ezdu1afve73/Hiiy8et3WyFm9FvKXgl6A5ePAgzzzzDABbt25lzpw5dO/e/YTz1N4R8ZZaPRIUffr0oVKlSowaNQqA7du3c/vttyvkRUKQZvwSkH379jFixAhSUlKoVasWU6ZMoWPHjtovXySEacYvZ+zYLH7+/PnEx8fz8ssvM3jwYDIzM0lMTPS2OBE5rXKZ8ZvZLcBoIAtIc859UB7XlbLh9/vZuHEjb7zxBvXr12fOnDlcf/31hce1eCsS2s56xm9mM8xsh5llFhtPNLMvzSzHzIYXDDvgR6ASkHv25YqXnHOkpKQAMHfuXMaMGUNGRsZxoQ9avBUJdYG0elKB4/5eb2YxwCSgIxAP3GVm8cCHzrmOwDDAH8A1xSODBw+mQoUK9OnTB4BDhw7x5JNPFj7FIyLh46yD3zm3Evi+2HALIMc5t8k5dwhIA7o45/IKju8GflbS7zOz/maWbmbpO3fuPNuyJMjy8vKYPHkyycnJXHDBBUyaNAnQN29FwlmwF3drAluKfM4FappZdzObCswCXirpB51z05xzCc65hOrVqwe5LDkbOTk5tGvXjgEDBtCqVSsyMzN58MEHvS5LRAIU7MVdK2HMOefeBN4M8rWkjDz11FP8/Oc/54knnqBixYokJydz3333YZb/X68Wb0XCW7CDPxe4vMjnWGBbkK8hZWjjxo2MHj0agE6dOjFlyhRq1qx53Dlq74iEt2C3etYAdcystplVBHoA84J8DSkDBw8exO/307hxYwBee+015s2bd0Loi0j4C+RxztnAauAaM8s1s77OuSPAQ8BiYCPwhnMuKzilSllZtWoVsbGx+Hw+Dh06BEDPnj2pUKGCZvciEeisWz3OubtOMr4QWHjWFUm52bt3L8OHD2fKlClcccUVLFiwgF//+teYmbZcEIlg2qsnCvl8Pho2bMhDDz3E9u3bGTJkCKNGjaJy5cpelyYi5UDBH2W2bt2K35//HbpGjRrx9ttv07x58+PO0VM7IpFNm7RFCeccSUlJxMfHAzB+/HjWrFlzQuiDntoRiXQK/iiwadMmrrzySvr378++ffsAGDZsGBUrVlTIi0QhBX8EO3r0KC+88AINGjRg165dTJ06lby8/N0ztOWCSPRS8Ecgn8/Hxo0buemmmxgyZAht27YlOzub/v37F377VkSil4I/whw+fLjwi1hffvklr776KvPnzyc2NrbwHC3eikQ3BX8ESU9PL9wbv0uXLmRnZ9OzZ88TZvlq74hENwV/BPjhhx+4/vrrad68OZ999hkAf/3rX/nFL36hkBeREyj4w9zcuXOpV68ea9asYcCAAezZswfQ4q2InJyCPwz5fD62bNlC165d6datG9WqVePvf/87L730ElWrVvW6PBEJcQr+MHP06FH8fj/x8fEsWbKE8ePHk56eTsuWLQvP0eKtiJyKgj+MfP7554UB37p1a7Kyshg6dCjnnnvuceepvSMip6LgDwMHDx7k5ptvplGjRqSnpwOwaNEirrzySoW8iJwxBX+I+/jjj2natCkrV67k3nvv5d///jegxVsROXsK/hDk8/nYv38/jzzyCDfccAP79u1jwYIFzJw5k2rVqnldnoiEOW3LHIL8fj+zZs1i06ZN/OEPf2DcuHFceOGFhce1eCsigVDwh5Dvv/+eYcOGAVChQgVWrFhBmzZtTjhP7R0RCYRaPSHAOUfXrl2pVq0a06dPByAnJ4ebb75ZIS8iQafg91hWVha33HILb7/9Nq1atWL9+vWAFm9FpOwo+D1wbPF2+PDhNG7cmMzMTJKSkvjoo49o2LCh1+WJSIRTj98Dfr+f1NRUvvnmG+677z7Gjx9P9erVC49r8VZEylK5BL+ZdQX+B6gBTHLOLSmP64aarVu38tBDDwFQuXJlVq5cyU033XTCeWrviEhZOutWj5nNMLMdZpZZbDzRzL40sxwzGw7gnJvrnOsH9AbuDKjiMJSXl0enTp2IjY1l7ty5QH5vv02bNgp5ESl3gfT4U4HEogNmFgNMAjoC8cBdZhZf5JQnC45HjX/84x+0a9eOBQsW0K5dO3JycgAt3oqId846+J1zK4Hviw23AHKcc5ucc4eANKCL5RsPvOucW3v25YYHn8/H4cOHGTduHA0bNmTdunVMnz6dZcuWcdVVV3ldnohEuWD3+GsCW4p8zgWuBwYC7YGqZna1c25K8R80s/5Af4ArrrgiyGWVL7/fz7x58/jss8/o3r07L730EpdddlnhcS3eioiXgh38VsKYc879H/B/p/pB59w0YBpAQkKCC3Jd5eLHH39k9OjRAHz77bfMmTOH7t27n3Ce2jsi4qVgP8efC1xe5HMssC3I1wg5zjnuuOMOqlSpwrPPPgvAd999x+23366QF5GQE+zgXwPUMbPaZlYR6AHMC/I1QsqGDRto27Ytc+bMoUmTJqxatQrQ4q2IhK5AHuecDawGrjGzXDPr65w7AjwELAY2Am8457KCU2ro8Pl87Nmzh4cffpgmTZqwYcMGpkyZwpo1a7jhhhu8Lk9E5JTMudBrpyckJLhjb5oKNXl5ecTExFC9enV27drF/fffz+jRo4/bJ9/n82mmLyLlzswynHMJpztPe/WcgfXr19OqVSsA6tatS0ZGBpMnTz7h5SgKfREJZQr+Uvjpp59o3bo1jRs35tNPPwVg1apVNGnSRCEvImFHwX8a77//Pg0bNmTVqlXcd9997Nq1C9DirYiELwV/CXw+H7t376Zv3760a9cO5xzLly9nxowZXHzxxV6XJyISEAV/Mc45/H4/9erV45VXXmHYsGFs2LCBdu3aFZ6jb96KSDjTfvxFfPHFFzz22GMAxMbGsmjRIho3bnzCeWrviEg404yf/H3ymzZtSr169ViwYAEAGRkZWrwVkYgUlcF/LMz37NnDiBEjuPrqq8nMzGTw4MHs2LED0OKtiESuqGz1+P1+qlSpwtixY9mzZw89e/Zk1KhR1K5d2+vSRETKXFTN+J1zzJo1C4DHHnuM66+/nrVr1zJr1qzjQl+LtyISyaIm+B977DEqVKjA//7v/xaOLVq0qPBViEWpvSMikSwqWj1vvvkmqampVKpUifHjxzNo0CBCcY8iEZHyENEz/r1799K7d29uv/124uLiWLt2LQ8//LDXZYmIeCoig9/n87FixQoaNWrErFmzGDlyJKtXr6ZevXqAevgiEt0iblvmAwcOcN5552FmXHXVVcyaNYuWLVsGuUIRkdATldsyb926lebNmwNw//33s27dOoW+iEgxERP8Pp+P2NhYMjMzAZgyZQqVK1fWEzoiIsVEXKsHwMz01I6IRJ2obPWIiMjpRWTw66kdEZGTi8jgV19fROTkIjL4RUTk5BT8IiJRRsEvIhJlFPwiIlFGwS8iEmVC8gtcZrYT+MbrOk7hEuDfXhdxCqovMKovMKovMIHUV8s5V/10J4Vk8Ic6M0svzbfjvKL6AqP6AqP6AlMe9anVIyISZRT8IiJRRsF/dqZ5XcBpqL7AqL7AqL7AlHl96vGLiEQZzfhFRKKMgv8MmNlmM9tgZuvM7OxfGBBEZjbDzHaYWWaRsYvNbKkqxNHAAAADmUlEQVSZfVXwz5+HWH0+M9tacB/XmdmvPartcjN738w2mlmWmQ0qGA+J+3eK+kLl/lUys0/NbH1Bff6C8dpm9knB/XvdzCqGWH2pZvbPIvevsRf1Fakzxsw+M7N3Cj6X+f1T8J+5ts65xiH0OFgqkFhsbDiw3DlXB1he8NkrqZxYH8DzBfexsXNuYTnXdMwR4FHnXD2gJTDAzOIJnft3svogNO7fQaCdc64R0BhINLOWwPiC+uoAu4G+IVYfwONF7t86j+o7ZhCwscjnMr9/Cv4w55xbCXxfbLgL8ErBn18BupZrUUWcpL6Q4Jz71jm3tuDPP5D/L19NQuT+naK+kODy/Vjw8dyC/zigHfC3gnEv79/J6gsZZhYL/A8wveCzUQ73T8F/ZhywxMwyzKy/18WcwqXOuW8hPzyAGh7XU5KHzOzzglaQZ62oY8wsDmgCfEII3r9i9UGI3L+CNsU6YAewFPga2OOcO1JwSi4e/p9V8fqcc8fu39iC+/e8mf3Mq/qAF4ChQF7B52qUw/1T8J+ZG51zTYGO5P+1u43XBYWpl4GryP/r97fAc14WY2aVgTnAYOfcPi9rKUkJ9YXM/XPOHXXONQZigRZAvZJOK9+qily4WH1mVh8YAVwLNAcuBoZ5UZuZdQJ2OOcyig6XcGrQ75+C/ww457YV/HMH8Bb5/0MPRdvN7DKAgn/u8Lie4zjnthf8C5kHJOHhfTSzc8kP1decc28WDIfM/SupvlC6f8c45/YAH5C/FnGRmZ1TcCgW2OZVXccUqS+xoIXmnHMHgRS8u383Ap3NbDOQRn6L5wXK4f4p+EvJzC4wsyrH/gzcBmSe+qc8Mw/oVfDnXsDbHtZygmOhWqAbHt3Hgn5qMrDROff/ihwKift3svpC6P5VN7OLCv58HtCe/HWI94E7Ck7z8v6VVN8XRf5P3cjvn3ty/5xzI5xzsc65OKAH8J5zriflcP/0Ba5SMrMryZ/lA5wD/MU5N9bDkgAws9nALeTv6LcdeBqYC7wBXAH8C/itc86TBdaT1HcL+W0KB2wG7j/WUy/n2loDHwIb+G+P9Y/k99E9v3+nqO8uQuP+NSR/8TGG/EnkG865UQX/rqSR30b5DLinYHYdKvW9B1Qnv62yDnigyCKwJ8zsFuAx51yn8rh/Cn4RkSijVo+ISJRR8IuIRBkFv4hIlFHwi4hEGQW/iEiUUfCLiEQZBb+ISJRR8IuIRJn/D1BGZuiyDwAIAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Nmax = 41;\n", "N = np.arange(2,Nmax);\n", "c = np.zeros(Nmax-2);\n", "for n in N:\n", " x = np.linspace(-1,1,n);\n", " V = np.vander(x,increasing=\"True\");\n", " c[n-2] = np.linalg.cond(V);\n", "\n", "pl.semilogy(N,c,'k-+');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A better way to go about interpolating with polynomials is via Lagrange interpolation. Define the Lagrange polynomial $L_i(x)$ to be $1$ when $x=x_i$ and is zero at all the other nodes, i.e.,\n", "$$L_i\\bkt{x_j} = \\delta_{ij}$$\n", "We then have\n", "$$p_n(x) = \\sum_{i=0}^n f_i L_i(x)$$\n", "Since we want $L_i(x)$ to vanish at all $x_j$, where $j \\neq i$, we have $L_i(x) = c_i \\prod_{j \\neq i}\\bkt{x-x_j}$. Further, since $L_i\\bkt{x_i} = 1$, we get that $c = \\dfrac1{\\prod_{j \\neq i} \\bkt{x_i-x_j}}$. Hence, we see that\n", "$$L_i\\bkt{x} = \\prod_{j \\neq i} \\bkt{\\dfrac{x-x_j}{x_i-x_j}}$$\n", "If we call $l_i(x) = \\prod_{j \\neq i} \\bkt{x-x_j}$ and $w_i = \\prod_{j \\neq i} \\bkt{\\dfrac1{x_i-x_j}}$, we see that $L_i(x) = w_i l_i(x)$.\n", "\n", "Further, if we set $l(x) = \\prod_{j=0}^n \\bkt{x-x_j}$, we see that $l_i(x) = \\dfrac{l(x)}{x-x_i}$, and hence $L_i(x) = \\dfrac{w_il(x)}{x-x_i}$ and hence we see that\n", "\\begin{align}\n", "p_n(x) = l(x) \\bkt{\\sum_{i=0}^n \\dfrac{w_i f_i}{x-x_i}} \\,\\,\\, \\spade\n", "\\end{align}\n", "Note that $\\spade$ is an attractive way to compute the Lagrange interpolant. It requires $\\mathcal{O}\\bkt{n^2}$ work to calculate $w_i$'s, followed by $\\mathcal{O}(n)$ work to compute the interpolant for each $x$. This is called as the ***first form of Barycentric interpolation***.\n", "\n", "What about updating when a new interpolation node $x_{n+1}$ is added? There are only two steps involved.\n", "\n", "- For $i \\in\\{0,1,\\ldots,n\\}$, divide each $w_i$ by $\\bkt{x_i-x_{n+1}}$. Cost is $n+1$ flops.\n", "- Compute $w_{n+1}$ for another $n+1$ flops.\n", "\n", "Hence, we see that the Lagrange interpolant can also be updated at $\\mathcal{O}(n)$ flops.\n", "\n", "The above barycentric formula can be made even more elegant in practice. Note that the function $1$ gets interpolated by any polynomial exactly. Hence, we see that\n", "$$1 = l(x) \\bkt{\\sum_{i=0}^n \\dfrac{w_i}{x-x_i}}$$\n", "which gives us that\n", "$$l(x) = \\dfrac1{\\bkt{\\displaystyle\\sum_{i=0}^n \\dfrac{w_i}{x-x_i}}}$$\n", "Hence, we obtain the ***second form of Barycentric interpolation***\n", "\\begin{align}\n", "p_n(x) = \\dfrac{\\bkt{\\displaystyle\\sum_{i=0}^n \\dfrac{w_i f_i}{x-x_i}}}{\\bkt{\\displaystyle\\sum_{i=0}^n \\dfrac{w_i}{x-x_i}}} \\,\\,\\, \\club\n", "\\end{align}\n", "Again note that it is fairly easy to incorporate a new intepolation node at a cost of $\\mathcal{O}(n)$. Below is the Lagrange interpolation using the original form." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c492f8186032456da9768f64229216f6", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=25, description='nnodes', max=45, min=5, step=2), Output()), _dom_classe…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def function(x):\n", "# f = np.abs(x)+x/2-x**2;\n", " f = 1.0/(1+25*x*x);\n", "# f = np.abs(x+0.3) + np.abs(x-0.2) + + np.abs(x*x*x*x-0.8);\n", " return f;\n", "\n", "def Lagrange(xnodes,x,i):\n", " f = 1;\n", " nnodes = np.size(xnodes);\n", " for j in range(0,i):\n", " f = f*(x-xnodes[j])/(xnodes[i]-xnodes[j]);\n", " for j in range(i+1,nnodes):\n", " f = f*(x-xnodes[j])/(xnodes[i]-xnodes[j]);\n", " return f;\n", " \n", "def Chebyshev(nnodes,xplot):\n", " # Chebyshev node interpolation\n", " xnodes = np.cos(np.arange(0,nnodes)*np.pi/(nnodes-1));\n", " fnodes = function(xnodes);\n", " fplot = 0;\n", " for i in range(0,nnodes):\n", " fplot = fplot + fnodes[i]*Lagrange(xnodes,xplot,i);\n", " return xnodes, fnodes, fplot;\n", "\n", "def Uniform(nnodes,xplot):\n", " # Uniform node interpolation\n", " xnodes = np.linspace(-1,1,nnodes);\n", " fnodes = function(xnodes);\n", " fplot = 0;\n", " for i in range(0,nnodes):\n", " fplot = fplot + fnodes[i]*Lagrange(xnodes,xplot,i);\n", " return xnodes, fnodes, fplot;\n", "\n", "def Bernstein(n,xplot):\n", " xnodes = np.linspace(-1,1,nnodes);\n", " fnodes = function(xnodes);\n", " fplot = 0;\n", " for i in range(0,nnodes):\n", " fplot = fplot + sp.stats.binom.pmf(i,nnodes-1,0.5+0.5*xplot)*function(xnodes[i]);\n", " return fplot;\n", " \n", "nplot = 1001;\n", "xplot = np.linspace(-1,1,nplot);\n", "f_actual = function(xplot);\n", "@interact\n", "def inter(nnodes=(5,45,2)):\n", " xnodes, fnodes, fplot = Chebyshev(nnodes,xplot);\n", "# xnodes, fnodes, fplot = Uniform(nnodes,xplot);\n", " # fplot = Bernstein(nnodes,xplot);\n", "\n", " error = f_actual-fplot;\n", " print(np.amax(np.abs(error)))\n", " pl.plot(xplot,f_actual,'-');\n", " pl.plot(xplot,fplot,'r');\n", " pl.rcParams[\"figure.figsize\"] = [16,4];" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }